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Introduction 


Our  study  of  helium  clusters  was  motivated  by  the  desire  to  understand  the  scaling  of  the  unusual 
properties  of  bulk  4He,  a  quantum  liquid,  in  finite  size  systems  as  one  goes  from  the  macroscopic 
regime  to  the  regime  of  molecular  dimensions.  This  is  fully  in  the  spirit  of  general  cluster  research, 
namely  to  develop  our  understanding  of  how  the  transition  from  molecular  to  bulk  systems  (or  vice 
versa)  is  reflected  in  the  internal  structure  and  dynamics  of  finite  size  aggregates.  The  unique  feature 
of  helium  is  its  dominant  quantum  behavior,  resulting  from  a  low  mass  and  weak  interatomic  binding 
energy.  Clusters  of  helium  are  therefore  very  weakly  bound  van  der  Waals  species,  whose  properties 
were  expected  to  be  dominated  by  zero  point  delocalization  effects.  During  this  grant  period,  \vc 
devoted  our  attention  exclusively  to  clusters  of  4He,  which  are  bose  systems.  These  are  more  strongly 
bound  than  the  fermionic  species  3HeN  ,  and  are  also  easier  and  cheaper  to  study  experimentally. 
Furthermore,  analogy  with  the  bulk  behavior  suggests  that  any  superfluid  effects  if  present  will  occur 
at  considerably  higher  and  therefore  experimentally  more  readily  accessible  temperatures  for  the  4  He 
species.  In  addition  to  the  helium  clusters,  we  also  applied  our  Monte  Carlo  techniques  to  clusters  of 
molecular  hydrogen,  which  for  J=0  are  also  Bose  clusters.  These  are  more  strongly  bound  than 
clusters  of  helium,  yet  are  still  very  delocalized  by  chemical  standards  and  offer  an  intriguing 
possibility  of  a  new  superfluid. 

Goals  of  original  research  plan 

1)  To  understand  the  size-dependent  scaling  of  superfluid  behavior  or  analogous  collective  effects  in 
clusters  of  4HeN.  As  a  preliminary  step  this  involved  analysis  of  Bose-Einstein  condensation  in  a 
weakly  interacting  bose  cluster. 

2)  Despite  much  phenomenological  progress  in  the  understanding  of  superfluidity  in  bulk  helium  II.  a 
molecular  description  for  the  characteristic  excitations  of  buik.h^lium  found  only  in  the  superfluid  state 
was  still  missing.  By  developing  a  truly  microscopic  theory  of  collective  excitations  in  these  quantum 
clusters  involving  accurate  ground  and  excited  state  wavefunction  information,  we  aimed  to  achieve 
new  insight  into  the  atomic  dynamics  underlying  the  superfluid  state  in  bulk  He  II  by  identifying  and 
analyzing  the  behavior  in  finite  sized  clusters  . 

3)  Determination  of  feasible  experimental  probes  of  the  cluster  dynamics,  in  particular  of  charged 
species  or  of  non-dissociative  molecular  probes.  This  goal  is  further  related  to  the  more  long  term  aim 
of  employing  the  unusual  physical  properties  of  these  quantum  clusters  to  modify  and  control  the 
course  of  chemical  reactions  of  embedded/attached  species  at  ultra-low  temperatures. 

4)  Development  of  Monte  Carlo  methods  to  provide  accurate  ground  and  excited  state  wavefunctions 
for  the  4HeN  clusters.  This  original  aim  was  extended  to  deal  also  with  the  more  strongly  bound  (H.)s 
species,  which  have  more  solid-like  character. 
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Goals  achieved  during  grant  period 


Our  central  achievements  attained  during  the  grant  period  are  the  following: 

1)  Establishment  of  a  new  quantum  liquid  drop  theory  for  the  collective  excitations  of  these  Bose 
clusters,  (papers  1, 2, 4).  Together  with  the  accurate  ground  state  wavefunctions  described  below, 
this  led  to  calculations  of  the  compressional  excitation  spectrum  for  L  =  0  (spherically  symmetric)  and 
L  =  1  (dipole)  symmetry,  for  a  range  of  cluster  sizes.  The  size  scaling  of  the  excitation  spectrum 
showed  the  onset  of  a  roton  minimum  at  sizes  N  ~  100,  leading  to  the  important  physical  conclusion 
that  cluster  of  size  N  >  100,  corresponding  to  diameters  R  >  10  A  can  support  superfluid  flow. 

2)  Development  of  accurate  Monte  Carlo  methods  for  ground  state  wavefunctions  of  general  quantum 
clusters  (papers  5,  7,  8,  11).  These  consisted  of  both  variational  and  (exact)  diffusion  Monte  Carlo 
techniques.  While  the  basic  'unguided'  Metropolis  sampling  of  variational  cluster  wavefunctions  had 
been  previously  employed  by  the  nuclear  phys:  :s  community,  we  improved  the  accuracy  and  sampling 
techniques  considerably,  by  developing  new  variational  wavefunction  forms  and  using  guiding 
functions  to  optimize  sampling  at  small  interparticle  separations.  This  resulted  in  an  unprecedented 
precision  of  -  5%  in  density  profiles  in  the  interior  of  the  cluster,  and  led  to  an  unexpected  discovery 
of  a  large  collinear  contribution  to  the  structure  of  the  He3  trimer  (paper  7).  Application  of  diffusion 
Monte  Carlo  to  these  weakly  bound  atomic  systems  is  new,  and  was  performed  selectively  to  calibrate 
the  variational  results  (paper  8).  The  diffusion  Monte  Carlo  algorithm  was  modified  to  employ  a 
constant  ensemble  size  (paper  11),  and  also  to  allow  modifications  of  the  quantum  force  in  regions  of 
small  inter-particle  dis*  nee  (paper  8).  High  accuracy  for  relatively  large  cluster  sizes  was  achieved  as 
a  result  of  refining  the  original  diffusion  Monte  Carlo  algorithm  (papers  8,  11).  These  exact 
calculations  show  structure  in  both  density  profiles  and  pair  distribution  functions  which  is  compatible 
with  some  very  weak  hard-core  packing  effects.  No  such  structure  is  observed  at  the  variational  level. 

3)  Development  of  a  quantum  theory  for  atomic  and  molecular  impurities  in  helium  clusters  (papers  7. 

1 1).  This  has  so  far  been  restricted  to  ground  state  impurities,  and  consists  of  a  variational  approach  to 
the  new  cluster  containing  the  impurity,  in  which  the  latter  interacts  pairwise  with  the  helium  atoms. 

The  variational  wavefunction  is  extended  to  include  pairwise  correlation  terms  between  the  impurity 
and  the  helium. 

This  approach  was  applied  first  to  a  H2  impurity  attached  to  HeN,  N  =  2  - 19,  and  structural 
analysis  of  the  resulting  mixed  cluster  ground  state  made  (paper  7).  The  lighter  impurity  H2  is 
extensively  delocalized  throughout  the  cluster,  with  a  peak  in  the  vicinity  of  the  diffuse  surface  region  ■ 
Quantitative  analysis  showed  that  the  H2  is  however  still  largely  in  the  interior  of  the  cluster  for  this 
range  of  sizes.  The  He13  species  is  unique  in  the  extent  to  which  it  expels  the  H2,  suggesting  an 
unusual  structural  stability  which  may  be  associated  with  an  ieosahedral  unit.  This  is  particularly 
significant  in  light  of  the  absence  of  any  energetic  magic  numbers  for  neutral  helium  clusters. 
Subsequent  unpublished  calculations  for  a  D2  impurity  show  quite  a  different  situation,  with  the 
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heavier  D2  species  considerably  more  localized  at  the  cluster  center,  and  with  no  peak  in  or  near  the 
surface  region.  The  helium  distribution  is  complementary  to  the  D2  distribution,  showing  a  peak  away 
from  the  cluster  center.  However  for  the  trimer  species  D2He2  the  central  peak  corresponding  to  near 
collinear  configurations  still  persists. 

Both  variational  and  diffusion  Monte  Carlo  approaches  were  then  applied  to  the  analysis  of  SF6  in 
HeN  (paper  11).  This  is  a  much  more  complex  system,  due  to  both  a  greater  binding  of  the  impurity 
with  He,  and  also  to  the  marked  anisotropy  of  the  He-SF6  interaction  potential.  New  trial 
wavefunctions  were  developed  which  faithfully  incorporated  the  main  features  of  the  potential 
anisotropy.  Calculations  were  made  for  clusters  SF6  with  N  varying  from  1  (the  'dimer')  to  N=499. 
The  ground  state  structures  of  these  systems  showed  pronounced  localization  of  the  helium  around  a 
centrally  located  SF6  impurity.  This  localization  occurs  in  sequential  shells,  with  the  first,  nearest 
neighbor  shell  containing  about  22-23  atoms  at  a  density  comparable  to  that  of  bulk  solid  He  at  -100 
atmospheres.  The  shells  further  away  from  the  impurity  contain  successively  larger  numbers  of  He 
atoms,  and  for  sizes  up  to  N=499  are  still  at  densities  lower  than  that  of  bulk  liquid  helium.  However 
it  is  interesting  that  the  exchange  energy  for  substitution  of  one  He  by  the  impurity  (paper  7)  appears 
nevertheless  to  have  saturated  by  N=  111,  despite  that  fact  that  the  structure  has  not  converged  to  what 
we  expect  for  an  SF6  impurity  in  bulk  helium.  For  SFg  we  also  analyzed  the  spectral  shifts  of  the  o~ 
vibrational  absorption  lines  due  to  the  instantaneous  dipole-induced  dipole  (IDID)  interactions  with  the 
surrounding  'solvent'  helium  species.  This  was  motivated  by  recent  experimental  measurements  of 
such  shifts  for  SF6  in  HeN  (Goyal,  S.,  Schutt,  D.  L.,  and  Scoles  G.,  Phys.  Rev.  Lett.  69,  933 
(1992))  and  in  (H2)N,  (Goyal,  S.,  Schutt,  D.  L.,  Scoles  G.,  and  Robinson,  G.  N„  Chem.  Phys.  Lett 
196,  123  (1992).  We  find  a  red-shifted  absorption,  which  increases  with  size  to  a  value  of  -  0.93 
cm'1  for  N=1 1 1,  with  a  half-width  of  -  0.25  cm'1 .  This  is  somewhat  smaller  than  the  experimental 
value  (~  1.5  cm'1)  which  also  differs  from  our  calculation  by  being  split  into  two  components.  The 
most  likely  reason  for  this  difference  is  that  in  the  experiment  the  clusters  are  not  in  the  ground  state, 
but  gain  a  considerable  amount  of  angular  momentum  from  pick-up  of  the  impurity.  Large  amounts  of 
angular  momentum  cause  considerable  centrifugal  distortions,  as  we  summarize  in  4)  below,  and  may 
cause  the  SF6  to  be  located  in  an  asymmetric  position  which  can  give  rise  to  a  line  splitting,  e.g.,  at  a 
cluster  surface.  This  possibility  can  be  investigated  with  Monte  Carlo  techniques  based  on  trial 
wavefunctions  combining  the  features  of  these  impurity  functions,  and  of  the  excited  rotational  states 
described  below.  Such  studies  are  planned  for  future  work. 

4)  Development  of  quantum  theoretical  approaches  to  excited  states  for  the  collective  modes  (papers 
2,  5,  10).  This  began  with  a  variational  approach  to  excited  compressionai  states  which  was  based  on 
the  Feynman  operator  approach  (papers  2,  5).  The  first  four  excited  compressionai  states  for  L  =  0 
were  calculated  variationally  for  N  =  240,  maintaining  orthogonality  to  lower  states  by  a  generalized 
Gramm-Schmidt  procedure.  These  results  showed  a  significant  lowering  of  the  compressionai 
energies  relative  to  both  classical  estimates  based  on  the  conventional  macroscopic  liquid  drop  model, 
and  also  to  our  new  quantum  liquid  drop  model.  We  have  recently  extended  the  variational  approach 
to  excited  states  of  overall  rotation  of  the  cluster,  employing  a  diffeicnt  approach  from  the  Feynman 


3 


excitation  operator  (paper  10).  New  trial  functions  which  are  eigenfunctions  of  total  angular 
momentum  and  which  are  physically  motivated  to  have  rotational  rather  than  vibrational  character, 
were  developed,  and  employed  in  both  variational  and  fixed  node  diffusion  Monte  Carlo  studies  of 
He7  and  of  its  more  strongly  bound  counterpart  (H2)7.  The  rotational  energies  are  considerably  lower 
than  the  compressional  energies  studied  previously  (papers  2,  5).  The  rotationally  excited  states  for 
He7  were  found  to  be  metastable  with  respect  to  dissociation  for  L  >  2,  while  (H,)7  is  metastable  only 
for  L  >  6.  Both  species  showed  very  large  oblate  centrifugal  distortion  at  low  L  values,  which 
developed  to  diffuse  toroidal  distributions  at  large  L.  The  distribution  of  the  cluster  surface  clearly 
changes,  which  has  possible  consequences  on  the  distribution  and  spectral  shifts  of  impurities  such  as 
SF6  (see  3)  above). 

5)  The  ground  and  excited  state  Monte  Carlo  techniques  developed  for  helium  clusters  have  also  been 
applied  to  clusters  of  H2  (J  =  0),  which  is  also  a  Bose  system  (papers  6,  10,  12).  We  introduced  an 
additional  element  of  employing  'shadow  wavefunctions’  for  fictitious  particles  representing  lattice 
sites  here  in  order  to  allow  for  more  rigid  structures.  The  primary  goal  of  this  extension  of  our  helium 
studies  is  to  determine  whether  a  liquid-like  ground  state  exists  for  (H2)N  for  N  small,  and  if  so, 
whether  these  clusters  display  similar  superfluid  behavior  to  HeN.  Our  preliminary  results  (paper  6) 
showed  that  the  smallest  clusters  (N  <  7)  are  extensively  delocalized.  More  recent  results  based  on  a 
more  accurate  trial  wavefunction  and  on  subsequent  diffusion  Monte  Carlo  calculations  (paper  1 2) 
show  that  clusters  up  to  N  =  33  still  show  strong  delocalization,  although  weak  features  characteristic 
of  solid  like  close  packed  structures  are  now  also  apparent.  We  have  analysed  the  structure  of  both 
these  and  the  HeN  clusters  in  the  body  fixed  frame  by  computing  principal  moments  of  inertia,  thereby 
avoiding  the  orientational  averaging  implicit  in  usual  Monte  Carlo  sampling  for  finite  systems. 
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Superfluidity  in  Helium  Clusters 


M.V.  Rama  Krishna' 1  and  K.B.  Whaley* 


20.1  Introduction 

Theoretical  and  experimental  investigations  of  clusters  has  been  an  active  and  growing  area 
of  research  for  the  past  several  years.  Yet,  much  of  our  current  understanding  of  the  structure, 
dynamics,  and  energetics  of  clusters  is  based  on  essentially  “classical”  clusters.  By  this  we 
mean  that  although  quantum  mechanics  is  important  in  the  description  of  their  electronic 
structure,  it  does  not  play  a  role  in  the  statistical  behavior  of  the  atoms  (or  rather,  of  the  nuclei) 
themselves.  However,  when  the  atoms  in  the  cluster  are  light,  such  as  H  or  He,  quantum 
mechanics  plays  a  significantly  different  role.  It  is  thus  important  to  study  quantum  clusters 
in  order  to  understand  what  role  quantum  statistical  effects  may  play  in  clusters,  and  how  the 
uniquely  quantum  phenomena  such  as  superconductivity  and  superfluidity  are  modified  in 
finite  systems.  With  this  in  mind  we  will  discuss  clusters  of  4He,  which  obey  Bose  statistics. 
The  goals  are  threefold:  I)  to  understand  how  superfluidity  manifests  itself  in  helium  clusters, 
2)  to  determine  the  energy  level  spectra  of  these  clusters,  and  3)  to  establish  experimental 
probes  of  these  clusters. 


20.2  Transition  Temperatures 

Bulk  liquid 4 He  is  known  to  undergo  a  phase  transition  from  a  normal  to  a  superfluid  state 
at  about  2.17  K.  This  phase  transition  is  characterized  by  a  nearly  logarithmic  divergence  of 
the  heat  capacity,  and  by  the  fact  that  the  superfluid  phase  can  flow  through  fine  capillaries 
with  zero  viscosityJ1*  Quantum  statistics  is  the  key  to  this  effect.  For  example,  liquid  3He, 
which  is  made  of  fermions,  exhibits  superfluid  behavior  only  at  a  much  lower  temperature  of 
about  1  x  10~3  K.  This  difference  is  not  simply  a  mass  effect,  but  results  from  the  need  of 
the 3 He  nuclei  to  pair  to  form  effective  Bose  particles.  Other  symmetry  related  properties  are 
significant  as  well.  For  example,  if  one  models  liquid  4 He  as  a  non-interacting  Bose  gas  one 
finds  that  this  model  exhibits  a  cusp  in  the  heat  capacity  curve  at  about  the  same  temperature 
where  the  experimental  heat  capacity  curve  exhibits  divergency.^ 

Given  these  observations,  and  their  relation  to  the  quantum  behavior  of  the  particles  mak¬ 
ing  up  liquid  4He,  we  wish  to  understand  how  the  phase  transition  is  modified  in  4Heyv  clusters 
due  to  finite-size  effects.  Although  there  is  no  longer  a  true  phase  transition  in  clusters,  it  is 
reasonable  to  use  the  same  ideas  that  have  been  previously  used  to  understand  the  phase  tran- 


*  Department  of  Chemistry,  University  of  California,  Berkeley,  CA  94720. 

•Present  address:  Box  955  Havemeyer  Hall,  Columbia  University,  New  York,  NY  10027-6948. 


258 


M.\  .  Rama  Krishna  and  h'.li.  Whaley  /  Superfluidity  in  Hr  Clusters 


sition  in  bulk  liquid  4Hc.  Consequently,  we  begin  by  discussing  the  “transition  temperatures” 
and  the  condensate  fractions  tor  clusters. 

The  Bose-Einstein  (BE)  condensation  temperature  of  a  non-interacting  Bose  gas  of density 
P  is  given  by*3* 


=  4.015  x  10-’V/3  K.  (20.1) 

For  liquid  helium,  Eq.  (20. 1 )  gives  I'h  /,•  =  3. 1 3  K,  which  is  0.96  K  larger  than  the  experimental 
value.11,  -1  We  used  Eq.  (20.1)  to  calculate  I'm-:  of  helium  clusters,  with  the  modification  that 
we  subtract  0.96  K  from  the  computed  values  so  that  for  sufficiently  large  clusters  we  recover 
the  experimental  bulk  value  of  7\  correctly.  These  corrected  Tlt  arc  given  in  Table  20.1.  The 
densities  of  the  clusters  arc  determined  using  p  =  .3/(4  th-3),  where  the  r0  are  the  calculated 
unit  radii.141 

One  can  use  Ginzburg-Landau-Pitacvskii  (GLP)  theory  to  estimate  the  transition  temper¬ 
atures  of  the  interacting  system.151  This  is  a  phenomenological  theory  in  which  the  free  en¬ 
ergy  density  is  expanded  in  terms  of  the  order  parameter,  which  here  is  an  effective  complex 
wavefunction  of  the  fluid.  The  expansion  is  valid  only  when  the  order  parameter  is  small 
and  the  coherence  length,  which  is  the  length  scale  over  which  the  order  parameter  changes, 
is  large.  Consequently,  this  theory  is  applicable  only  when  the  temperature  T  of  the  fluid  is 
dose  to  T\.  In  the  original  mean-field  GLP  theory,  the  expansion  coefficients  were  functions 
of  integer  powers  of  (7\  -  T).151  Such  a  mean-field  approach  neglects  fluctuations,  which  are 
very  important  close  to  T\ .  The  modem  version  of  this  theory  due  to  Mamaladze  employs  a 
modified  free  energy  density  that  accounts  for  the  fluctuations  of  the  order  parameter  near  T\ 
by  taking  the  temperature  dependence  of  the  expansion  coefficients  from  experiments.16, 71 
The  successful  applications  of  this  modified  GLP  theory  to  the  prediction  of  transition  tem¬ 
peratures  of  helium  films  and  pores17, 81  gives  us  confidence  regarding  its  utility  in  the  case 
of  clusters.  Clearly  this  theory  can  not  predict  critical  exponents  as  these  are  put  in  via  the 
expansion  coefficients.  However,  these  exponents  can  also  be  obtained  from  first  principles 
starting  from  the  GLP  theory  and  using  the  renormalization  group  theory  of  Wilson.  For  our 
purposes  here  this  is  not  necessary. 

When  applied  to  a  spherical  cluster,  the  modified  GLP  theory  predicts  the  transition  tem¬ 
perature  7\  to  be  given  by1?1 

rA  =  7\6-^  K.  (20.2) 

II  i 

where  R  is  the  radius  of  the  cluster  in  X,  T*  is  the  transition  temperature  of  bulk  liquid  helium, 
and  T\  is  that  of  the  cluster.  Of  course,  although  the  transition  temperature  is  a  sharply  defined 
quantity  in  macroscopic  systems,  in  clusters  we  expect  that  there  is  rather  a  temperature  range 
over  which  the  transition  to  the  superfluid  state  takes  place.  The  temperatures  predicted  above 
essentially  give  the  location  of  the  peak  values  of  the  rounded  heat  capacity  curves  of  the 
clusters. 

To  calculate  the  condensate  fraction  in  these  clusters  at  T  =  0  K,  let  us  use  the  model  of 
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Table  20.1:  Transition  temperatures  and  zero-temperature  condensate  fraction  in  helium  clus¬ 
ters.  (R  is  given  in  X  ,  The  and  7\  in  Kelvin.) 


N 

R 

I 'be 

T\ 

9cC 

20 

7.4 

1.1 

0.9 

32.7 

40 

8.8 

1.4 

1.2 

27.1 

70 

10.2 

1.6 

1.4 

22.6 

112 

11.8 

1.6 

1.5 

21.2 

240 

14.7 

1.8 

1.7 

17.1 

728 

20.9 

1.9 

1.9 

15.0 

10000 

47.8 

2.1 

2.1 

12.2 

Liq.  He 

x: 

2.17 

2.17 

9.2 

an  imperfect  Bose  gas  for  the  clusters.  This  is  given  bv^ 


Ao 

.V 


3\A 


«3/y/2 


(20.3) 


=  1 -6.148p1/2.  (20.4) 

where  a  =  2.556  X  is  the  experimentaJ  scattering  length  of  the  helium  atoms.  The  percent 
condensation  %C  is  simply  Ao/A'  x  100.  These  results  are  given  in  Table  20.1. 

We  see  from  Table  20.1  that  the  theoretical  estimates  of  The  and  T\  agree  remarkably 
well  even  for  a  cluster  as  small  as  twenty  atoms,  and  that  the  bulk  transition  temperature  is 
depressed  by  only  about  0.5  K  in  Hc24o-  We  also  find  that  the  condensate  fraction  approaches 
that  of  the  bulk  fluid  rather  rapidly.  Note  that  the  condensate  fraction  is  decreasing  as  T\  is 
increasing.  This  indicates  that  the  strong  interactions  between  particles  in  the  denser  (larger) 
clusters  are  depleting  the  zero  temperature  condensate,  while  increasing  the  transition  tem¬ 
peratures. 


20.3  Collective  Excitations 

Another  quantity  related  to  superfluidity  is  the  excitation  spectrum  J25  Fora  microscopic  un¬ 
derstanding  of  superfluidity,  it  is  important  to  understand  how  the  collective  excitation  spec¬ 
tra  of  clusters  change  as  a  function  of  cluster  size,  and  how  they  approach  that  of  the  bulk 
fluid.  With  this  goal  in  mind,  consider  helium  clusters  as  quantum  liquid  drops  of  radius  R 
and  uniform  interior  density  po.  The  density  wave  excitations  of  the  droplet  are  characterized 
by  the  quantum  numbers  (/.  w)  and  n.  The  momenta  kin  of  these  excitations  arc  given  by  the 
boundary  condition  ji(kinR )  =  0.  Within  this  liquid  drop  model  a  harmonic  analysis  of  the 
collective  vibrations  of  the  cluster  gives^ 10, 1 1 ' 

fi2k2 

‘n  2mS,(k,n)ui„ 


(20.5) 
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Excitation  Spectra 


Figure  20.1:  Excitation  spectra  of  /  =  0  and  1  modes  obtained  using  Eq.  (20.5)  for  N  =  20 
(circle),  70  (square),  and  240  (triangle)  clusters.  The  tilled  symbols  give  1=0  spectra  and  the 

open  ones  give  l  =  1  spectra.  The  solid  line  ( - )  is  the  Bijl-Feynman  excitation  spectrum 

of  bulk  He  II. 


where  .5)  is  the  structure  function  of  the  cluster,  defined  as  the  Fourier  Bessel  transform  of 
the  density-fluctuation-density-fiuctuation  correlation  function,  and  uln  is  the  normalizauon 
factor  for  the  spherical  Besse!  functions.  Complete  details  of  the  theory  and  calculations  are 
given  in  Refs.  10, 1 1. 

Equation  (20.5)  is  the  finite  cluster  analog  of  the  Bijl-Feynman  excitation  spectrum  for 
bulk  liquid  helium J1*  It  represents  the  compressional  vibrational  excitation  energies  of  the 
cluster,  and  in  the  bulk  limit  corresponds  to  the  phonon  spectrum  of  liquid  He.  To  get  a  picture 
of  these  spectra  one  needs  to  compute  the  structure  functions  St.  Monte  Carlo  random  walk 
simulations  for  (  =  0  and  I,  and  N  =  20, 70,  and  240  were  performed.  The  spectra,  together 
with  the  Bijl-Feynman  excitation  spectrum  of  liquid  helium,  are  shown  in  Fig.  20.1.  We  see 
that  the  spectrum  of  the  clusters  evolves  toward  that  of  the  bulk  fluid  rather  rapidly.  The 
pronounced  dip  at  k  «  2  X-1  in  the  liquid  helium  spectrum  is  known  as  the  roton  region. 
The  He?o  cluster  already  shows  such  a  roton  structure  at  about  2  X~‘  and  the  spectrum  of 
He24o  strongly  resembles  that  of  liquid  helium.  The  validity  of  these  results  are  confirmed 
by  a  more  general  theory  based  on  Bijl-Feynman  wavefunctions  for  the  excited  states  of  the 
clusters.^1" 

Since  the  excited  states  of  a  many-body  system  play  an  important  role  in  both  the  ther¬ 
modynamics  and  dynamics  of  the  system,  the  strong  resemblance  of  the  excitation  spectra  of 
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He/v  cluster.  10  the  bulk  fluid  is  indicative  that  these  clusters  will  also  exhibit  similar  ther¬ 
modynamic  and  dynamic  behavior.  It  is  reasonable  therefore  to  expect  helium  clusters  of 
about  tOO  atoms  to  undergo  a  normal  superfluid  transition  strongly  resembling  that  of 
liquid  helium.  There  are  additional  arguments  for  making  this  connection  between  the  excita¬ 
tion  spectrum  and  superfluidity.  For  example,  Bogoliubov  first  showed  that  the  phonon-roton 
spectrum  of  liquid  helium  is  a  result  of  both  interactions  between  Bose  particles  and  the  pres¬ 
ence  of  the  Bose-Einstein  condensate.* 121  (See  also  Table  20.1  and  the  associated  remarks 
made  in  See.  20.2.)  For  the  simple  model  of  a  weakly  interacting  Bose  gas  he  obtained 


i  ( k )  = 


\frrV  r.v\ 


r  .v ' 
T~ } 


(20.6) 


where  ('  =  4ml/ tn  is  the  interaction  energy  assumed  to  be  constant  and  repulsive,  and  k  is 
the  linear  momentum  associated  with  'he  excitation.  This  yields  a  linear  part  (phonon  branch) 
at  small  momenta  and  a  quadratic  piece  (Ircc-particlc  branch)  at  higher  momenta,  and  thus 
reproduces  the  mam  features  of  the  phonon-roton  s|>ecirum  of  liquid  helium.  Furthermore, 
if  cither  ( '  or  /<  =  X/V  is  very  small,  then  the  spectrum  will  be  completely  frec-paruclc 
like.  Such  a  system  w ill  not  exhibit  superfluidity  since  the  minimum  velocity  needed  to  excite 
the  fluid,  known  as  the  Landau  critical  velocity,  vanishes  in  this  case.121  Other  liquids,  such 
as  water,  which  exhibit  collective  excitations*13'  at  large  k  values  (/■  >  1  X-!)  also  may 
not  exhibit  superfluidity  because  these  liquids  are  stable  only  at  high  temperatures.  In  this 
regime  the  thermal  excitations  dissipate  the  energy  of  the  moving  particles.  Although  the 
model  of  a  weakly  interacting  Bose  gas  is  not  quantitatively  appropriate  for  liquid  helium,  the 
essential  relationship  between  Bose-Einstein  condensation  and  the  phonon-roton  spectrum  is 
still  present.  Hence,  the  onset  of  a  phonon-roton  type  spectrum  fora  Bose  fluid  is  a  signature 
of  a  large  Bose-Einstein  condensate  and  of  superfluidity.  Based  on  this  argument  the  results 
presented  in  Fig.  20.1  give  evidence  that  clusters  of  about  70  atoms  should  be  superfluid  at 
sufficiently  low  temperatures. 

Recently,  path-integral  Monte  Carlo  simulations  have  been  used  to  compute  heat  capaci¬ 
ties  and  superfluid  densities  of  He*4  and  He^g  clusters  as  a  function  of  temperature.*14*  The 
peaks  of  the  computed  heat  capacity  curves  yield  “transition  temperatures”  of  1.6  and  1.8  K, 
respectively,  in  very  good  agreement  with  those  reported  in  Table  20.1  for  similar  sized  clus¬ 
ters.  The  path-integral  simulations  also  indicate  that  the  width  of  the  heal  capacity  maximum 
is  increasing  while  T\  is  decreasing,  with  decreasing  cluster  size.  Consequently,  we  antici¬ 
pate  the  phase  transition  to  be  completely  washed  out  in  clusters  of  about  20  atoms  or  less.  It 
will  be  interesting  to  pursue  the  study  of  superfluid  densities  as  a  function  of  temperature  for 
a  series  of  small  clusters,  to  see  if  some  of  these  clusters  arc  indeed  non-supcrfluid  even  at  0 
K. 


20.4  Detecting  Superfluidity 

Two  research  groups  have  been  making  pioneering  efforts  to  detect  superfluidity  in  free  helium 
clusters,  but  the  experimental  evidence  so  far  is  inconclusive.*13, 16*  This  stems  primarily 
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electron  helium  pseudopolenlial 


Figure  20.2:  New  cIcctron-He  pscudopoteniiai  of  Ref.  IS  ( - )  is  compared  with  the  best 

previous  potential.  Ref.  19  ( - ). 


from  the  difficulty  in  probing  these  extremely  weakly  bound  van  dcr  Waals  clusters,  which 
are  easily  dissociated  and  whose  internal  excitations  have  until  now  been  poorly  understood. 
Heat  capacity  measurements!’7'  on  bubbles  of  helium  confined  in  copper  foil  have  shown 
presence  of  the  superfluid  state  in  bubbles  of  radius  40-60  °A,  corresponding  to  N  >  104. 
These  experiments  provide  direct  evidence  for  the  depression  of  T\  and  rounding  of  the  heat 
capacity  peak  in  these  finite  systems,  although  quantitative  analysis  of  7*  is  complicated  by 
the  presence  of  the  confining  copper  matrix. 

The  possibility  of  binding  an  electron  to  the  surface  of  helium  clusters  was  considered  re¬ 
cently,  in  order  to  use  it  as  a  spectroscopic  probe  of  the  cluster.  A  simple  but  accurate  electron- 
He  interaction  potential  (pscudopoteniiai)  capable  of  reproducing  s-  and  p-wave  scattering 
phase  shifts  accurately  over  a  range  of  energies  is  needed  as  a  start.  Such  a  pscudopotential 
has  only  recently  been  developed.'18'  In  Fig.  20.2,  it  is  compared  with  the  commonly  used 
prior  pseudopotential  for  this  system.'19'  The  well  depth  of  about  300  K  is  almost  ten  times 
deeper  than  the  previous  pscudopotential.  The  phase  shifts  and  scattering  cross  sections  cal¬ 
culated  using  this  new  pscudopotential  are  given  in  Table  20.2.  These  reproduce  s-  and  p- 
wave  scattering  phase  shifts  to  within  1 .5%,  and  total  and  momentum-transfer  cross  sections 
to  within  3%  of  the  exact  values,'20'  over  a  range  of  electron  energies  from  0-16  cV.  This  is 
a  major  improvement  over  the  previous  pseudopotential,  which  gives  s-wave  phase  shifts  to 
within  only  38%,  and  yields  p-wave  phase  shifts  with  an  incorrect  sign. 

Using  this  very  accurate  pseudopotential  for  the  short-range  interaction,  and  a  polarization 
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Table  20.2:  Phase  shifts  and  cross  sccuons  calculated  using  the  new  pseudopoiential.  The  % 
errors  given  in  the  parentheses  arc  calculated  using  the  estimated  values  given  in  Tables  III 
and  IX  of  Ref.  20.  The  total  and  momentum-transfer  cross  sections  are  obtained  from  the  s- 
and  p-wave  phase  shifts  calculated  using  the  c-Hc  \  pseudopoiential  of  Ref.  18  and  Eqs.  (7) 
and  ( 1 0M 1 2)  of  Ref.  20. 


E(eV) 

Phase  Shifts  (radians) 
s-wave  p-wavc 

Cross  Sections  (x  10  16  cm') 
Total  Momentum  transfer 

0.136 

-0.12700  (-0.9) 

0.00311  <  0.9) 

5.656  (-1.8) 

5.929  (-1.8) 

0.544 

0.26581  (  0.1) 

0.01318  (0.5) 

6.118  (  0.2) 

6.695  (  0.3) 

1.224 

-0.40785  (  1.4) 

0.03069  (  0.2) 

6.265  (  2.7) 

7.107  (  2.5  :> 

2.177 

-0.54667  (1.5) 

0.05526  (  0.1) 

6.151  (  2.6) 

7.1.67  (  2.3) 

3.401 

-0.67824  (  1.5) 

0.08572  (-0.4) 

5.860  (  2.3) 

6.898  {  2.0) 

4.898 

0.80070  {  1.0) 

0.120.35  (-0.4) 

5.473  (  1.3) 

6.417(1.0)  | 

6.666 

0.91352(0.8) 

0.15719  (-1.0) 

6.047  (  0.7) 

5.802  (  0.4) 

8.707 

-1.01691  (  0.1 

0.19439  (-0.8) 

4.616  (  -0.0) 

.6.131  (0.6) 

11.020 

-1.11137  (  -0.4) 

0.23043  (  -0.0) 

4.200  (-0.4) 

4.460  (  0.2) 

13.605 

-1.19750  (-0.7) 

0.26422  (  0.6) 

3.808  (-0.3) 

3.830  (  0.0) 

16.462 

1.27604  (-0.7) 

0.29508  (  0.6) 

3.445  (-0.1) 

3.26)  (0.1) 

potential  for  the  long-range,  we  obtain  a  complete  c-Hc  v  potential.118*  The  energy  levels  of 
the  excess  electron  may  now  be  determined  for  clusters  of  various  sizes.  One  of  the  interest¬ 
ing  results  of  this  calculation  is  that  it  takes  approximately  5  x  105  helium  atoms  to  barely 
bind  the  electron  with  a  binding  energy  of  about  0.04  cm-1.  Also,  the  electron  is  very  “dif¬ 
fusely”  bound  to  the  clusters  (see.  Fig.  20.3),  which  unfortunately  makes  it  insensitive  to  the 
internal  structure  and  dynamics  of  the  clusters.  When  the  cluster  is  large  enough  to  be  consid¬ 
ered  bulk  liquid  helium,  the  experimental  zero-field  energy  levels  of  the  excess  electron  are 
reproduced.* 1 8*  This  gives  us  confidence  in  the  accuracy  of  the  pseudopoiential  and  also  in  the 
calculated  energy  levels  of  the  excess  electron  on  clusters.  While  it  is  possible  in  principle  to 
determine  the  sizes  of  the  large  clusters  by  comparing  the  experimental  and  theoretical  values 
of  the  energy  of  the  photon  needed  to  barely  detach  the  excess  electron,  the  very  weak  binding 
does  not  make  this  a  very  useful  probe.  We  are  investigating  the  use  of  embedded  molecules 
as  alternative  indirect  spectroscopic  probes  of  the  cluster  size  and  dynamics.  Moreover,  the 
accurate  elcctron-He  \  pseudopoiential  now  available  makes  it  possible  to  model  e-He.v  scat¬ 
tering  experiments,  which  could  possibly  provide  insight  into  ways  to  excite  the  collective 
states  of  the  cluster,  and  thereby  yield  experimental  confirmation  of  the  calculated  excitation 
spectrum. 


20.5  Summary 

Theoretical  evidence  has  been  presented  for  the  existence  of  superfluid  helium  clusters  at  ex¬ 
perimentally  accessible  temperatures.  Large  clusters  of  several  hundred  atoms  appear  to  un¬ 
dergo  a  transition  to  a  superfluid  state  strongly  resembling  that  of  bulk  liquid  helium.  Exper- 
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Figure  20 .3:  The  square  of  the  wavefunction  of  an  excess  electron  attached  to  tne  surface  of 
various  He,v  clusters,  versus  the  distance  r  from  the  cluster.  The  cluster  size  N  is  indicated 
on  the  plots. 


imental  detection  of  such  a  superfiuid  state  remains  an  outstanding  problem.  Theoretically, 
it  remains  to  be  investigated  whether  small  clusters  with  N  <  20  will  be  superflutd  in  any 
sense  at  any  temperature,  and  what  atomic  motions  are  responsible  for  the  superfluid  dynam¬ 
ics.  Theoretical  efforts  are  also  necessary  to  determine  the  best  possible  experimental  means 
of  detecting  superfluid  behavior  in  these  extremely  weakly  bound  clusters. 
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Monte  Carlo  study  of  impurities  in  quantum  clusters:  H2  4HeN,  N=2-19 
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Variational  Monte  Carlo  techniques  are  employed  in  studying  4He  dusters,  with  and  without 
an  H,  impurity.  We  find  that  a  novel,  yet  simple,  analytic  nuclear  wave-function  form,  derived 
from  a  numerical  H2He  wave  function,  yields  high  accuracy  in  computed  ground-state  energies 
of  4HeN .  For  the  clusters  studied  here,  three  to  twenty  atoms,  energies  range  from  94%  to 
90%  of  the  exact  values.  Density  profiles  and  distributions  of  panicle  separation  are  also 
computed.  For  reasonable  computational  cost  (e.g.,  <  20  Cray/X-MP14  minutes  for  the 
largest  duster),  density  profiles  are  determined  for  the  first  time  to  high  statistical  accuracy  to 
within  0.5  A  or  less  of  the  duster  center.  The  density  profile  of  He,  is  found  to  possess  a 
uniquely  pronounced  peak  at  the  cluster  center  resulting  from  contributions  of  near-collinear 
atomic  arrangements.  We  also  study  the  effect  of  substituting  an  He  by  H2,  using  modified 
wave  functions  containing  products  of  pairwise  He-H,  terms.  For  all  duster  sizes  studied,  v.  * 
find  a  lowering  of  the  total  energy  upon  exchanging  an  He  for  an  H2.  The  exchange  energy 
increases  in  magnitude  with  increasing  duster  size,  yet  is  still  well  below  bulk  estimates  at 
N  ~  20.  Size  comparisons  with  the  pure  helium  dusters  show  very  little  change  upon  He/H2 
exchange,  e.g.,  the  rms  radii  differ  by  <2%  for  N  >  3.  Density  profiles  and  bond  distributions 
show  noticeable  differentiation  between  H2  and  He.  For  N>4,  the  peak  in  the  H2  density 
profile  is  not  at  the  cluster  but  does  remain  inside  the  duster.  This  peak  is  most  pronounced  for 
H.He,  j  implying  an  enhanced  resistance  to  H:  penetration  for  He,,. 

l 


I.  INTRODUCTION 

Clusters  of  rare  gases  constitute  the  simplest  van  der 
Waals  aggregates.  These  species  continue  to  provide  much 
stimulus  for  both  experimental  and  theoretical  investigation 
of  size-dependent  properties  in  the  regime  spanning  molecu¬ 
lar  and  bulk  systems. 1  Theoretical  study  is  especially  attrac¬ 
tive  because  the  interaction  potentials  are  so  well  known. 
The  most  weakly  bound  of  these  species  is  4HeN ,  in  which 
large  quantum  effects  are  apparent,  making  them  of  special 
fundamental  interest  for  the  understanding  of  size-depen- 
dent  behavior  of  quantum  systems.  Recent  predictions  for 
these  Bose  clusters,  based  upon  scaling  of  the  microscopic 
excitation  spectra2  and  the  calculation  of  the  phenomeno¬ 
logically  defined  superfluid  fraction,1  has  for  the  first  time 
provided  a  quantum  analog  to  recent  extensive  analysis  of 
phase  transitions  and  structural  features  of  classical  clus¬ 
ters.4 

While  these  weakly  bound  species  intrinsically  provide 
ideal  testing  grounds  for  new  theoretical  approaches  to  large 
quantum  aggregates,  contact  with  experimental  studies  has 
thus  far  proved  elusive.  The  ease  of  dissociation  and  the  li¬ 
quidlike  structure  of  helium  clusters,  which  facilitates  the 
absorption  of  foreign  species,  has  also  rendered  the  results  of 
scattering  experiments  ambiguous.5,6  Nevertheless,  there 
does  now  exist  a  considerable  body  of  experimental  data  on 
the  pickup  abilities  and  ionization  patterns.6  It  has  become 
clear  that  rather  than  attempting  a  direct  probe  of  the  dy¬ 
namics  and  excitations  of  these  quantum  clusters,  indirect 
probes  via  spectroscopic  studies  of  embedded  species,  prefer¬ 
ably  without  the  additional  complications  introduced  by  dis¬ 
sociation  channels,  will  provide  more  useful  information.  In¬ 
frared  spectroscopy  of  small  molecules  in  argon  clusters  has 


yielded  information  on  the  location  of  the  foreign  species.' 
The  same  information  is,  in  principle,  available  for  the  quan¬ 
tum  clusters  of  helium,  and  furthermore,  analysis  of  spectro¬ 
scopic  line  shapes  should  give  information  on  the  coupling  to 
the  internal  cluster  excitations  and  whether  these  are  collec¬ 
tive  or  single-partide-iike. 

In  this  paper  we  undertake  a  fully  quantum-mechanical 
study  of  an  embedded  molecule  in  helium  clusters.  We  em¬ 
ploy  variational  Monte  Carlo  ( VMC)  methods  in  studying 
the  structural  and  energetic  effects  of  adding  a  foreign  spe¬ 
cies,  using  a  new  wave-function  form  based  on  accurate  pair 
potentials.  The  techniques  are  applied  here  to  H2  in  He* 

( from  here  on  4HeN  is  assumed ) ,  a  choice  motivated  by  both 
the  availability  of  high-quality  He-He  and  H2-He  poten¬ 
tials,  and  also  by  prior  experimental  observations  of  H2  in 
bulk  helium.*  We  calculate  the  ground-state  energy  of  clus¬ 
ters  with  N<20,  with  and  without  an  H2  species  attached, 
and  thereby  the  He/H2  exchange  energy.  Analysis  of  the 
density  profiles  resulting  from  the  optimized  wave  functions 
shows  that  the  embedded  Hj  is  extensively  delocalized 
throughout  the  cluster,  with  a  small  peak  in  the  concentra¬ 
tion  beneath  the  surface.  Such  structural  analysis  clearly  de¬ 
scribes  the  location  of  the  foreign  species. 

Very  few  species  are  even  metastable  in  bulk  helium  be¬ 
cause  of  its  closed-shell  configuration  and  inert  nature. 
However,  there  has  been  continued  interest  in  the  analysis  of 
impurities  in  bulk  helium  because  of  the  properties  of  these 
both  as  nucleation  centers,9  and  as  a  source  of  information 
on  the  induced  response  and  spatial  structure  of  the  bulk 
medium.  The  only  theoretical  work  to  date  for  species  other 
than  the  isotopic  impurity  5He  has  employed  paired  phonon 
analysis  in  the  hypemetted-chain  approximation  with  Len- 
nard-Jones  potentials.10  As  pointed  out  by  Kfirten  and  Ris- 
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tig  in  Ref.  10,  these  techniques  are  subject  to  large  errors. 
However,  it  has  proven  difficult  to  proceed  beyond  this  level 
of  description,  either  by  VMC  or  Green's  function  Monte 
Carlo  ( GFMC),  because  of  the  relatively  small  contribution 
of  the  impurity  to  the  total  energy  in  the  bulk  system.  Analy¬ 
sis  of  the  cluster  energies  with  and  without  impurities  as  a 
function  of  cluster  size  provides  a  means  of  approaching  the 
bulk  impurity  problem  by  extrapolation  and  therefore 
yields  a  new  microscopic  approach  to  the  impurity  probe 
problem  for  bulk  helium  as  well  as  for  clusters. 

The  energies  computed  here  for  the  pure  HeN  clusters 
compare  favorably  with  previous  VMC  and  obtain  a  high 
percentage  of  the  exact  diffusion  Mcnte  Carlo  ( DMC )  ener¬ 
gies.  In  our  wave  functions  only  two-body  correlation  fac¬ 
tors  are  used.  One-body  factors,  often  employed  in  cluster 
calculations,  are  not  found  to  be  necessary  here,  and  three- 
bodv  factors  are  seen  to  become  significant  only  for  N  S  20. 
Details  of  the  physical  motivation  of  our  wave-function  form 
and  a  thorough  description  of  the  various  Monte  Carlo  ap¬ 
proaches  employed  are  presented  elsewhere. 11  The  remain¬ 
der  of  this  paper  is  organized  as  follows.  Section  II  contains  a 
description  of  our  theoretical  approach,  with  a  brief  sum¬ 
mary  of  the  guided  and  unguided  Metropolis  walks  used 
here,  in  addition  to  description  of  the  potentials,  wave-func¬ 
tion  forms,  and  optimization.  Results  for  pure  and  mixed 
clusters  are  presented  and  discussed  in  Sec.  III.  Conclusions, 
together  with  a  prognosis  for  further  studies  of  excitations 
and  dynamics  of  clusters  with  foreign  species,  are  presented 
in  Sec.  IV. 

It.  THEORETICAL  APPROACH 

The  quantum  ground  states  are  studied  by  seeking  accu¬ 
rate  manv-body  nuclear  wave  functions.  That  is.  in  place  of 
unknown  eigenfunctions  of  the  nuclear  Schrodinger  equa¬ 
tion.  we  obtain  wave  functions  for  which  the  energy 
E—  <  VP!// ivk)/('I<l'I'>  is  minimized.  For  the  systems  we 
treat  here,  the  Hamiltonian  H  is,  in  atomic  units, 

-£(2/n,)-‘V?  +  2FV  (!) 

'  KJ 

where  a  pairwise  potential  is  considered  sufficient  for  weakly 
interacting  clusters.12 

For  the  pure  He  clusters  m,  =  mHt  and  Vtj  =  V  is  the 
accurate  and  widely  employed  interatomic  potential  deter¬ 
mined  by  Aziz.  McCourt,  and  Wong  in  1987.13  This  poten¬ 
tial  has  a  slightly  deeper  well  than  the  previous  one  deter¬ 
mined  in  1979, 14  and  yields  a  bound  state  for  the  He-He 
dimer  at  —  10' 3  K.  For  clusters  containing  H2,  vibration  of 
the  foreign  molecule  is  neglected.  This  should  be  quite  rea¬ 
sonable  for  these  weakly  bound  clusters  since  the  vibrational 
dependence  of  the  H;-He  potential  is  significant  only  at 
small  separations  R  ( R  is  the  distance  from  the  H;  center  of 
mass  to  He ) .  In  addition,  the  H2  -He  potential  is  very  nearly 
isotropic.13  Expanding  FH..He(/?,0)  in  the  Legendre  poly¬ 
nomials  P0  (cos  6)  and  Fj  ( cos  6)  ( 9  is  the  angle  between  the 
i  n  termolecuiar  axis  and  R )  yields  an  anisotropic  component 
V,(R)  with  a  well  depth  an  order  of  magnitude  smaller  than 
that  of  yn(R).  Therefore,  as  a  further  yet  still  accurate  ap¬ 
proximation,  we  choose  =  K,.  Thus,  in  the  Hamilto¬ 


nian  we  employ,  H:  acts  as  a  single  parude  with  twice  the  H 
atom  mass  and  interacts  with  He  via  K0. 

The  potential  is  a  Lennard-Jones  plus  van  der  Waais 
fit  to  the  ab  initio  data  of  Ref.  1 5. 

f4e[(cr/R)'2  -  (a/R)”],  R<R., 

V(R\  = 

0  —  CtR  —  CSR  —  C.0R  ~ l0.  R>R.. 

(2) 

where  Rr  is  chosen  such  that  the  two  forms  in  Eq.  1 2 )  take 
essentially  the  same  value  at  this  point.  The  van  der  Waais 
parameters  are  taken  from  Ref.  15.  Table  I  lists  the  param¬ 
eters  describing  V0  and  Fig.  1  compares  with  the  ab  inno 
data  points  (the  root-mean-square  deviation  is  1  9*^  ) .  and 
with  the  He-He  potential. 

A.  Wave-function  form  and  optimization 

Once  a  wave-function  form  is  chosen,  parameters  are 
varied  to  minimize  the  energy,  E.  Here,  V  takes  a  transla- 
nonaily  invariant  product  form. 

T( Rl  =  </  i  r  ),  3) 

-•  / 

where  for  N  particles.  R  is  a  3N’-dimenstonal  vector  specify¬ 
ing  the  particle  locations,  and  r ,  =  t;  —  r,  |.  The  function 
<!>„ ,  analogous  to  the  /2  factor  defined  in  previous  studies  of 
atomic  c lusters. 21  vanes  with  differing  particle  pairs,  i  e  . 
He-He  or  H,  -He.  Note  that  a  one-body  factor /,  (r, )  is  not 
employed.  Such  a  factor  introduces  undesired  center-of- 
mass  motion.  In  this  event,  the  associated  translational  ener¬ 
gy  must  be  subtracted  ( unless  N  is  large),  adding  to  compu¬ 
tational  cost.  Alternatively,  translational  invanance  may  be 
maintained  by  replacing  r,  by  r,  —  Rcm  ,  where  Rcn,  is  the 
center-of-mass  vector,  yielding  a  sum  of  interparticle  separa¬ 
tions.  In  this  case,  therefore./,  is  actually  a  many-bodv  fac¬ 
tor.  A  common  motivation  for  including/,  is  its  use  in  bind 
ing  the  cluster  and  in  describing  the  more  diffuse  regions 
Since  our  two-body  factor  is  derived  from  an  exact  dimer 
( H,  He)  wave  function,  in  which  particular  attention  is  paid 
to  the  long-range,  diffuse  tail,  part  of  the  wave  function, :  1  we 
prefer  not  to  modify  this  by  addition  of  a  one-body  retaining 
factor.  If  the  cluster  is  bound,  this  should  be  reproduced  by  a 
pairwise  wave  function,  provided  it  is  of  sufficiently  high 
accuracy.  Our  new  approach  differs  from  the  original  one  of 
Pandharipande,  Pieper,  and  Wiringa,14  which  constrained 
the  asymptotic  form  of  tl>,t  by  the  relation  between  the  pair 
distribution  function  g(  r)  and  the  speed  of  sound  in  the  bulk 
system,  thereby  necessitating  a  one-body  factor  to  yield  a 
bound  finite  system.  Therefore,  although  a  one-body  factor 


TABLE  I.  Parameters  of  the  H.-He  interaction  potential.  P,.  Eq  i  Z  ■ 


Parameter 

Value 

Umu 

107 

4  189  712 

hartrees  i  h  > 

a 

5  671  145 

Bohr 

R, 

8.351  10 

Bohr 

C. 

4018 

H  (Bohr)* 

c. 

55.69 

h  (Bohr)" 

C,  o 

I  031.0 

h  (Bohr)"' 

-i  n  ojfndii  v  o  *<rai9/  .mouruies  n  hosiers 


FIG  1  H. -He  and  He-He  poientials  The  solid  line  jhow«  the  He-He  po- 
lennal  oi  Ref.  I J.  and  the  dotted  Itne  shows  the  H.  -He  potential.  K„ .  Eq. 

1 2 1  3nd  Table  I.  The  ab  initio  data  points  for  K, ,  given  in  Ref.  1 5.  are  indi¬ 
cated  by  the  sol  d  circles. 

has  proved  useful  in  earlier  studies,  it  is  extraneous  for  our 
purposes.  We  also  do  not  employ  a  three-body  factor  here 
because  of  the  high-quality  wave  functions  we  obtain  at  the 
two-body  level  of  wave-function  complexity.  We  shall  refer 
to  this  point  again  in  Sec.  III.  The  factors  are 

ii'„ ( r )  ss  /'explP(u)  +ar“],  (4) 

« 

Pi  u)  =  y  aku\  u  —  r  •  '.  (5) 

The  ij  subscripts  are  omitted  from  the  parameters  for  conve¬ 
nience:  however,  parameters  describing  different  interac¬ 
tions.  i  e..  V  and  Kn,  are  independent.  This  form  was  de¬ 
duced  from  the  numerical  solution  to  the  Schrodinger 
equation  employing  the  H-  -He  potential  given  in  Eq.  i  2 ) . " 

The  wave-function  form  may  now  be  considered  as  the 
product  of  a  long-range  and  a  short-range  factor.  At  large 
interparticle  separations  where  the  polynomial  P  is  roughly 
constant,  the  behavior  of  is  determined  by  a,  b,  and  a.  The 
long-range  factor,  /expfa^),  is  very  similar  to  that  em¬ 
ployed  in  previous  work,2'16""  with  the  exception  that  a  is 
now  a  variable  parameter  which  when  optimized  for  our 
wave  functions  is  different  from  unity.  Note  that  for 
b  =  -  ( N  -  1 ) ' '  and  a  =  1,  our  form  reduces  to  the  one- 
parameter  model  asymptotic  form  for  a  two-fragment  bound 
system.  ■*  The  remaining  factor.  exp(/*(u)  ],  is  most  impor¬ 
tant  at  small  separations  and  can  be  considered  as  a  general¬ 
ization  of  the  McMillan  form  used  for  bulk  He.20  Given  the 
highly  repulsive  behavior  of  the  potential  at  small  particle 
separations,  an  increase  in  wave-function  complexity  for  de¬ 
scribing  the  short-range  interactions  should  be  especially 
useful.  We  therefore  include  all  integer  powers  up  to  five  in 
P.  The  study  of  a  further  improved  two-body  form,  which 
yields  an  entirely  new  description  of  the  behavior  at  small 
particle  separations  tailored  to  the  functional  form  of  the 


repulsive  pan  of  the  potential,  will  be  considered  in  another 
work." 

The  major  motivation  in  employing  such  wave  func¬ 
tions  is  our  belief  that  considerable  improvement  over  pre¬ 
vious  two-body  forms  can  be  made,  to  the  extent  that  one  can 
efficiently  obtain  high  accuracy  without  the  use  of  a  complex 
three-body  factor.  The  general  questions  of  what  accuracy  is 
achievable  at  the  two-body  level  of  theory  and  the  extent  to 
which  our  wave  functions  obtain  this  accuracy  are  addressed 
in  a  subsequent  publication. ' 1  A  second  motivation  for  pro¬ 
ducing  such  accurate  wave  functions  is  the  need  for  high 
accuracy  when  considering  a  weakly  bound  impurity  spe¬ 
cies. 

Below  we  outline  our  wave-function  optimization  ap¬ 
proach  Full  details  are  given  in  Ref.  1 1.  We  optimize  the 
parameters  by  minimizing  fluctuations  of  the  local  energy. 
EL='V~  'ATI',  about  a  reference  (or guess)  energy,  £*.  over 
a  fixed  set  of  points  in  the  3N-dimensional  space  R  The 
quantity  of  interest  is  then 


V 


V  [£t(R,)  -  £*  ]2  V(R  WofR,)!2 


£  I'KR.J/jf^R,)*2 


<*![//- f,]2!*) 

1 1  i  (  Q  ) 

(*!♦> 

wheic  the  second  equation  corresponds  to  an  infinite  num¬ 
ber  of  points  sampled  from  i'K,  i2.  The  essence  of  this  "fixed- 
sample"  approach21  is  that  multidimensional  integrals  are 
approximated  as  summations  over  points  in  a  distribution 
corresponding  to  an  initial  (unoptimized)  wave  function. 
4V  Parameters  are  converged  to  optimum  values  (in  a  local 
minimum  sense)  with  a  conjugate  gradient  technique." 

The  reference  energy  E„  determines  the  emphasis  given 
to  minimizing  the  energy  vs  minimizing  the  variance.  If 
E„  4El,  minimizing  r  is  equivalent  to  minimizing  the  aver¬ 
age  local  energy,  where 


El 


£  £t(R1)|*(Rt>/¥o<*.>!J 

i  -  I _ _ _ 

£  |'P(R, )/'P0 (R, )l2 

<  »  l 


(7) 


However,  if  £*  =  £t,  then  minimizing  r  is  equivalent  to 
minimizing  the  variance  in  the  energy.  Correspondingly, 
choosing  intermediate  values  of  £„  mixes  energy  and  vari¬ 
ance  reduction.  As  we  are  primarily  interested  in  obtaining 
the  lowest-energy  wave  functions  possible,  we  generally 
choose  the  reference  energy  to  be  much  lower  than  EL 
Since  N  is  finite  (<5000  points),  the  summations  in  Eq 
(6)  may  be  poor  approximations  to  the  desired  expectation 
value.  This  is  especially  true  when  'P  becomes  much  different 
than  Y0 ,  which  occurs  if  is  a  poor  initial  guess.  Therefore, 
we  have  found  it  useful  to  set  ♦<,  =  'P  and  to  generate  a  new 
ensemble  from  the  updated  4V  The  conjugate  gradient 
minimization  is  then  continued.  This  updating  procedure 
greatly  enhances  the  reliability  of  the  degree  of  energy  mini¬ 
mization  indicated  by  an  individual  fixed-sample  optimiz¬ 
ation  Generally,  one  or  two  updatings  have  proved  suffi- 
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cient  to  obtain  stable  and  converged  results.  The  optimized 
wave-function  parameters  are  given  in  the  Appendix. 

B.  Monts  Carlo  approach 

We  are  interested  in  computing  expectation  values  of  4», 
which  cannot  be  done  analytically.  Numerical  quadrature 
may  be  employed  but,  given  the  high  dimensionality  in¬ 
volved  here,  Monte  Carlo  integration  is  the  most  tractable 
approach.  Here,  we  employ  two  variants  of  the  Metropolis 
walk.23  The  first  is  the  commonly  used  “unguided”  walk  by 
which  points  are  sampled  from  i 't' j2.  Here,  unguided  refers 
to  the  fact  that  the  attempted  moves  underlying  the  walk  are 
completely  random.  Expectation  values  are  computed,  to 
within  statistical  error,  as 

A  =N-'  £  4(R,)~<'t'U  lSk>/<^l*l'>.  (8) 

i  -  I 

All  previous  VMC  cluster  studies  employed  unguided 
walks. i  The  second  approach  used  here  employs  a  walk 
guided  by  a  different  function,  'l/r  For  these  walks,  attempt¬ 
ed  moves  are  biased  towards  larger  values  of  the  sampled 
distribution#  (Walks  guided  by  ^  itself  yielded  poor 
convergence  at  small  particle  separations.11  )  Expectation 
values  with  respect  to  'V  now  require  a  weighting  procedure, 
namely, 

As  =  £  ^(R,)l'KR,)/'l'1!(R/)|2 

£  i^(R,)/^(R,)|2J  (9) 

Generally,  the  usefulness  of  a  guided  walk  arises  from  im¬ 
portance  sampling,  i.e.,  the  preferential  sampling  of  (impor¬ 
tant  )  regions  where  good  statistics  are  required  for  high  pre¬ 
dion  in  Here,  —  qp  1S  chosen  so  that  i^i2  is 
preferentially  sampled  in  re  cions  where  it  is  large.  The  guid¬ 
ing  function  is  also  chosen  decay  to  zero  more  slowly  than 
^  at  small  r  to  facilitate  convergence  in  this  domain,  where 
significant  contributions  to  the  energy  are  typical.  A  particu¬ 
lar  advantage  of  the  "guided/weighted"  walk  lies  in  the 
flexibility  of  choosing  which  regions  are  to  be  emphasized. 
For  example,  if  high  precision  is  desired  at  small  separations, 
choosing  >  1  in  this  domain  yields  the  desired  sam¬ 

pling  The  guided/weighted  walk  was  most  useful  for  the 
smaller  dusters  studied  here.  Detailed  comparisons  of  the 
various  walks  and  their  relative  merits  for  these  weakly 
bound  systems  are  presented  in  Ref.  1 1. 

As  a  final  point,  we  discuss  the  evaluation  of  statistical 
error.  The  primary  concern  is  that  averaged  auantities  be 
statistically  independent  so  that  the  computed  statistical  er¬ 
ror  is  unbiased.  In  our  approach  this  is  accomplished  by 
propagating  (ten)  independent  ensembles  of  (100)  points 
{R,}  yielding  (ten)  independent  Monte  Carlo  estimates. 
Ideally,  this  can  be  accomplished  by  a  random  selection  of 
points.  In  practice,  for  the  smallest  clusters,  the  ensembles 
are  not  initially  decorrelated  but  commence  with  different 
random  number  seeds.  The  ensembles  are  propagated  in  par¬ 
allel  by  a  Monte  Carlo  walk  until  decorrelation  between 
them  is  obtained.  Decorrelation  may  be  ascertained  from  the 
stability  of  the  statistical  error  in  the  average  of  the  (ten) 


Monte  Carlo  estimates.  Also,  if  the  distance  moved  in  con¬ 
figuration  space  by  each  member  of  all  the  ensembles  is  sev  - 
eral  times  larger  than  the  dimensionality  of  the  cluster,  then 
the  ensembles  can  be  assumed  to  be  decorrelated.  To  reduce 
this  initial  decorrelation  time,  ensembles  for  larger  clusters 
are  built  up  from  the  ensembles  of  the  next  smallest  cluster 
In  addition,  this  parallel  structure  is  also  maintained  when 
generating  .urge  ensembles  of  points  by  VMC  walks  dunng 
wave-function  optimization.  Given  that  each  ensemble  is 
different  from  the  others,  employing  these  ensembles  in 
Monte  Carlo  runs  yields  a  set  of  independent  results.  This 
enables  the  computation  of  an  unbiased  statistical  error 
This  method  of  evaluation  of  the  statistical  error  is  used  for 
all  expectation  values,  including  those  yielding  distribution 
functions  such  as  density  profiles. 

111.  RESULTS  AND  DISCUSSION 
A.  Pur*  H*m 

In  this  work  we  study  clusters  ranging  in  size  from  three 
to  twenty  panicles,  the  focus  being  the  ground-state  energy 
and  structure.  The  pnmary  concern  here  is  wave-function 
accuracy,  as  measured  by  the  energies  obtained.  Therefore, 
we  first  compare  our  pure  He  cluster  energies  with  those 
resulting  from  other  (sometimes  more  complex )  wave  func¬ 
tions  and  with  new  DMC  energies  which  are  exact. 

Table  II  compares  our  optimized  He  cluster  energies 
obtained  with  the  new  wave  function.  Eqs.  (3)-(5).  wuh 
those  resulting  from  other  work.  Exact  energies,  from  which 
percent  accuracies  are  obtained,  are  the  DMC  values  given 
in  Ref.  1 1 .  For  N  =  3-3,  we  see  that  our  computed  energies 
lie  below  previous  VMC  values  which  were  computed  em¬ 
ploying  both  one-  and  two-body  factors,  indicating  ihe  high 
quality  obtainable  with  the  our  two-body  wave  functions  for 
these  small  dusters.  Similar  accuracy  has  also  been  obtained 
with  two-body  wave  functions  for  N  =  3-7  24  For  He;o  our 
new  energy  reproduces  to  within  the  statistical  error  ( given 
in  parentheses)  our  previous  result  obtained  with  one-,  two-. 


TABLE  H.  He  cluster  per  particle  energies.* 


This  work 

Previous  VMC 

-  £/N  (K) 

%  of  Exact 

—  £/N  (K.i 

Ts  of  Exacr 

He, 

-0.0415(1) 

93.9 

-00388(9)* 

81  7 

He. 

-  0.1357(1) 

93.9 

-0132000)* 

9!  3 

-our 

96  0* 

He, 

-0.2505(2) 

93.3 

-  0.2440(  10)* 

911 

He, 

-  0  4838<  1 ) 

92.7 

He,, 

-  1.0545(6) 

He,. 

-  1  1290(7) 

90.3 

He„ 

-  1.510(2) 

89.3 

-  1  5)4(3)* 

897 

-  1  573(1)* 

96  7* 

*  Here  and  in  Table  HI,  the  statistical  error  representing  one  standard  del¬ 
ation  m  the  mean  is  shown  in  parentheses.  Exact  energies  for  the 
recent  potential  are  computed  by  DMC  in  Ref.  1 1 . 

"Reference  I*. 

Reference  16. 

’Since  Ref.  1 6  employs  the  previous  potential  of  Ref.  IS,  the  tract  ciurg,  i 
taken  to  be  the  GFMC  value  of  Ref.  25. 
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and  three-body  factors,  and  gives  89.5%  of  the  exact  DMC 
energy.  (This  is  somewhat  less  accurate  than  the  VMC  re¬ 
sult  of  Pandharipande,  Pieper,  and  Wirings,16  also  obtained 
with  one-,  two-,  and  three-body  factors,  for  the  earlier  poten¬ 
tial,  14  which  yielded  96.7%  of  the  exact  GFMC  energy.26 ) 
Given  the  simplicity  of  our  wave  function,  the  proximity  of 
our  energy  to  that  resulting  from  one  which  is  more  complex 
is  quite  satisfying.  However,  we  see  that  for  very  high  accu¬ 
racy,  three-  and  perhaps  higher-body  factors  are  required  for 
N  X  20.  In  this  case,  our  present  wave  functions  appear  to 
provide  a  two-body  component  of  high  quality. 

Overall,  Table  II  shows  that  the  accuracy  of  our  ener¬ 
gies,  vis-a-vis  those  of  DMC,  decreases  with  increasing  clus¬ 
ter  size.  Naturally,  this  is  expected  as  a  pure  “two-body” 
wave  function  should  be  less  appropriate  as  the  importance 
of  three-body  and  more  complex  interactions  increase,  as  is 
the  case  for  larger  clusters.  However,  it  is  pleasing  to  see  that, 
despite  the  constant  level  of  wave-function  complexity,  the 
(energy)  accuracy  does  not  degrade  rapidly.  As  the  cluster 
size  increases  from  7  to  20  atoms,  the  amount  of  the  DMC 
energy  obtained  decreases  by  only  3%.  In  addition,  the  total 
accuracy  remains  good  up  to  N  =  20,  with  90%  of  the  DMC 
energy  obtained  for  this  cluster.  Therefore,  the  wave-func¬ 
tion  form  employed  here  is  considered  to  be  of  sufficiently 
high  flexibility  for  the  study  of  the  H,HeN  clusters,  N<  19. 

Helium  cluster  density  profiles  are  presented  in  Fig.  2. 
The  data  were  obtained  by  binning  the  sampled  points  into 
bins  corresponding  to  AJ?  <0.03  A.  Note  that  the  bin  volume 
and  therefore  the  statistical  precisian  decrease  as  the  cluster 
center  is  approached,  as  indicated  by  the  fluctuations  in  p 
which  well  represent  statistical  error.  In  spite  of  this,  preci¬ 
sion  remains  high  up  to  small  distances  from  the  cluster  cen¬ 
ter.  For  He20 ,  for  which  uncertainty  in  p  is  greatest,  the  com¬ 
puted  statistical  error  in  each  bin  is  under  5%  from  0.5  to  12 
A.  Thus,  any  peculiar  behavior  at  the  cluster  center,  such  as 


FIG. 2.  Helium  density  profiles  for  pure  helium  cluster  tires  N  «•  3-5, 7, 13, 
20.  The  bin  size  u  <0.03  k. 


a  hole,  as  well  as  any  hard-core  structure  in  the  single-parti¬ 
cle  distribution  function,  is  mostly  precluded  by  our  data. 
The  computational  cost  of  the  He  calculations  presented 
here  is  very  reasonable,  ranging  from  4  Cray/X-MP14  min¬ 
utes  for  He,,  to  19  CPU  minutes  for  Hew.  These  profiles 
represent  the  first  Monte  Carlo  study  of  any  type,  VMC, 
GFMC,  or  DMC,  which  possess  this  high  level  of  statistical 
accuracy  in  the  interior  of  the  clusters.  As  reflected  in  Fig.  2, 
the  efficiency  of  computing  density  profiles  decreases  with 
increasing  cluster  size.  For  a  cluster  with  N  particles,  each 
point  sampled  yields  N  values  to  be  binned.  However,  the 
computational  cost  increases  as  N2,  causing  the  efficiency 
for  this  quantity  to  scale  as  N  * '. 

Figure  2  demonstrates  the  unique  character  of  the  He, 
density  profile.  Approaching  the  cluster  center,  the  He,  pro¬ 
file  is  similar  in  character  to  those  of  the  other  clusters  up  to 
the  shoulder  located  at  about  2.2  A.  After  this  point,  a 
marked  change  in  behavior  occurs  as  p  increases  only  slight¬ 
ly  from  about  2.2  A  to  1  A,  and  then  rises  rapidly  closer  in. 
These  characteristics  are  presented  in  greater  detail  in  Fig.  3. 
In  addition  to  the  density  profile,  contributions  to p  are  pre¬ 
sented  conditional  on  the  largest  angle  (&„„)  in  the  He, 
triangle  being  greater  than  120",  140",  and  160".  This  figure 
shows  that  the  rapid  rise  in  p  near  the  origin  is  due  solely  to 
near-collinear  arrangements  starting  at  ~  1 20*.  That 
density  near  the  center  arises  solely  from  near-collinear  ar¬ 
rangements  is  not  unreasonable  given  the  high  potential  en¬ 
ergy  for  atoms  near  the  cluster  center  (and  each  other)  in  a 
near-equilateral  configuration.  However,  the  rapid  increase 
in  p  caused  by  a  tendency  towards  collinearity  is  interesting 
given  the  large  difference  of  a  collinear  structure  from  that  of 


FIG.  3.  Conditional  density  profiles  of  He, .  The  dashed  line  is  the  density 
profile.  The  three  conditional  density  profiles,  solid  lines,  are  obtained  by 
binning  distances  from  the  cluster  center  only  when  the  atomic  arrange¬ 
ment  sampled  possesses  an  angle  greater  than  120*.  I«F.  and  IMF. 
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minimum  potential  energy,  an  equilateral  mangle. 

In  addition  to  the  density  profiles,  we  present  panicle 
separation  probability  density  functions  p(r)  (normalized 
such  that  Jp(r)dr=*  1)  in  Figs.  4  and  5-  These  plots  were 
obtained  by  binning  particle  separations  with  a  bin  size 
Ar<0.03  A.  The  curves  in  Fig.  4  show  a  trend  toward  de¬ 
creasing  diffuseness  as  the  cluster  size  in  .reases  from  three 
to  five  atoms,  consistent  with  increasing  density  and  also 
increasing  binding  energy  per  particle.  Upon  considering  the 
larger  clusters  ( Fig.  5 )  this  trend  reverses.  The  interparticie 
separations  increase  slightly  from  He,  to  He,  and  more  no¬ 
ticeably  upon  going  to  He,,  and  He,„.  The  HeN  root-mean- 
square  ( rms )  radii  presented  in  Table  IV  mirror  this  behav¬ 
ior,  which  has  been  previously  noted  by  Pandhanpande, 
Pieper.  and  Wirings. 16  What  is  most  interesting  and  novel 
here  is  the  appearance  of  structure  in  p(  r )  for  He, ,  and  He:o 
(Fig.  5 ) .  The  behavior  of  p  as  it  approaches  its  maximum 
clearly  differs  between  He,  and  He, , .  For  He,a  the  existence 
of  a  shoulder  at  a  particle  separation  smaller  than  that  corre¬ 
sponding  to  maximum  p  is  unambiguous.  However,  the  pres- 
enceofa  shoulder  mp  for  He,,  is  less  certain.  The  He. »  plots 
of  the  p  and  p  are  virtually  identical  to  those  of  He,,  (and 
are.  therefore,  omitted  for  clarity).  This  similarity  between 
He, ,  and  Heu  is  in  contrast  to  the  larger  differences  between 
HeN  andHeN  * ,  observed  for  smaller  N  (Fig.  4).  Therefore, 
given  the  slowly  changing  nature  of  p  at  N  a  1 3.  we  conclude 
that  the  onset  of  the  shoulder  in  p  is  quite  gradual.  Such 
behavior  does  appear  to  correspond  to  evolution  of  shell 
structure,  and  is  related  to  the  appearance  of  the  second- 
nearest-ncighbor  coordination  shell  in  the  pair  density  dis¬ 
tribution  function  n:  (0 ,r)."  Our  density  functions  p,  to¬ 
gether  with  the  absence  of  any  structure  in  the  single-particle 
distribution  function  (Fig.  2),  are  consistent  with  the  inter- 


r  (A) 


FIG  4  Particle  separatum  probability  deputy  function  pin  for  the  pure 
helium  clusten.  N  «  3-5.  The  function  pin  is  the  probability  density  of 
finding  sny  two  particles  separated  by  the  distance  r.  The  curves  ere  nor¬ 
malised  to  unity  end  bin  sizes  are  <0.03  A. 


FIG  5  Particle  separation  probability  density  function  p t  r»  for  the  pure 
helium  clusters.  N  =  '  13. 20.  The  normalization  and  bin  stzesare as  in  F:e 
a 


pretation  of  He  clusters  as  delocalized,  liquidities  ci  isters 
rather  than  as  crystalline-  or  moiecularlike  species. 

B.  Mixed  clusters  H,HeM 

Table  III  lists  the  cluster  energies  per  particle  ( H,  is 
considered  as  one  particle)  and  the  He/H,  exchange  energy 
as  a  function  of  size.  The  exchange  energy  ( referred  to  as  ihe 
chemical  potential  for  He/H,  exchange  in  Ref.  10)  for  an  V 
panicie  cluster  is  computed  as 

£„(N)  =  £(H,HeN_, )  -£(Hev).  50: 

For  the  H ,  He  van  der  Waais  complex,  numerical  solution  of 
the  nuclear  Schrddinger  equation  yields  the  ground  state 
bound  by  0.0246  K.  while  the  He  dimer  is  bound  by  onh 
— 10  K.'J  This  is  not  surprising  given  the  greater  *eil 
depth  of  the  H,  -He  interaction  relative  to  the  He-He  poten¬ 
tial,  cf.  Fig.  1.  As  a  consequence,  negative  exchange  energies 
are  observed.  The  exchange  energy  increases  monotonicaily 
with  increasing  cluster  size;  however,  our  value  at  N  =  20. 
—  1.68(4)  K,  is  still  only  a  fraction  of  Kiirten  and  Ristiz  s 
estimate  of  the  bulk  value.  —  20  K.10  Table  (II  also  lists  the 
unit  radii,  r0,  defined  by 

r0  =  ($<*,>),'J.V-’'\  ‘It. 

where  <£  2)  is  the  expectation  value  of  the  squared  distance 
from  the  cluster  center.  The  latter  is  defined  as 


i.e..  the  geometric  center.  Since  the  center  of  mass  is  bia,co 
towards  heavier  particles,  we  employ  the  geometnc  center  j> 
r  he  origin  in  comparing  structure  and  size  of  pure  and  mi  seJ 
clusters,  and.  most  importantly,  in  studying  the  relative  c- 
cation  of  H,  as  specified  by  the  density  profiles.  Unlike  tr.e 
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TA8LE  III  Per  particle  energies,  unit  radii,  and  He/H,  exchange  energy.* 


He„ 

Hj  Hr* . , 

N 

£/N  (K) 

'a  (A) 

£/N  (K.) 

r,  (A) 

£„  (N) (K) 

3 

-0  041  52(8) 

5  50 

-0  089  38(24) 

5.21 

-0.143  6(8) 

4 

-  0.135  7(2) 

4.42 

-0.200  5(1) 

4.51 

-0  259  2(9) 

5 

-  0  250  5(2) 

3.76 

-0.332  4(2) 

3.84 

-0.409  5(1!) 

? 

-  0483  8(1) 

3.45 

-0.575  1(5) 

3.43 

-0.639  1(36) 

13 

-  1.054  5(6) 

297 

-  U46  1(U> 

2.98 

-  1.191(16) 

14 

-  1  129  0(7) 

2.92 

-  1.224  2(4) 

2.96 

-  1.333(  1  1  ) 

:o 

-  1.5100(20) 

2.82 

-  1.594  3(8) 

2.86 

-  1.680(44) 

•Statistical  error  in  r„  is  less  than  0.003  A. 


per  particle  energies,  cluster  size  is  largely  unaffected  by 
He/H,  exchange.  As  seen  in  Table  III,  the  unit  radii  of  pure 
and  mixed  clusters  differ  by  only  4%  for  three  panicles  and 
by  2%  or  less  for  the  other  clusters. 

In  considering  the  structure  of  the  mixed  clusters,  in 
particular  the  degree  of  differentiation  between  H,  and  He. 
we  first  turn  t6  the  root-mean-square  distance  from  the  clus¬ 
ter  center.  flrm,  =  (( R J))1/J.  Values  for  He  and  H,,  ob¬ 
tained  by  integrating  over  the  density  profiles  shown  in  Figs. 
6  and  7,  are  listed  in  Table  IV.  These  density  profiles, 
p(R)  =p(R),  are  normalized  as  S4ttR  2p(R)dR  =  1,  and 
as  stated  above  the  cluster  center  is  chosen  as  the  origin.  The 
corresponding  values  for  the  pure  clusters  are  also  shown  for 
purposes  of  comparison.  As  is  the  case  for  the  pure  clusters, 
we  see  that  a  minimum  in  cluster  size,  as  measured  by  Rrmt, 
also  occurs  at  five  particles  for  the  mixed  clusters.  The  differ¬ 
entiation  between  H2  and  He  Rrm  values  within  a  given 
cluster  ts  somewhat  small  but  consistent.  For  ail  clusters 
studied  here.  H,  is  farther  from  the  cluster  center  on  aver¬ 
age.  with  values  of  Rm,  being  4.6%  to  11%  greater  than 
those  of  He. 

Differentiation  between  H,  and  He  is  also  evident  in  the 
normalized  density  profiles  presented  in  Figs.  6  and  7.  To 
quantify  this  differentiation  between  H2  and  He,  we  com¬ 
pute  the  probability  P  of  finding  H2  (or  He)  outside  the 
spherical  volume,  centered  at  R  =  0,  containing  half  of  all 
the  particles.  Values  of  P  for  both  species  are  also  listed  in 
Table  IV  These  values  indicate  that  the  degree  to  which  H, 
lies  "outside"  the  cluster  is  generally  not  large.  For  all  but 
the  14-particie  cluster,  60%  or  less  of  the  H2  density  lies 
outside  the  half-particle  dividing  radius.  The  exception, 
H.He,,.  for  which  66%  lies  outside  this  radius,  actually 
continues  a  trend  from  the  four-mrticle  mixed  cluster  which 
then  reverses  in  proceeding  to  H;  He,,.  This  suggests  a  spe¬ 
cial  structural  robustness  of  He,3  to  H2  penetration.  Pre¬ 
vious  OFMC  results  for  He  argue  against  the  presence  of 
"magic  numbers”  for  the  pure  clusters,  N<33.JS  However, 
the  effect  we  see  here  is  small  and  also  may  not  be  deducible 
from  the  energies  of  pure  clusters.  Whether  the  increase  in 
penetration  of  H2  seen  on  going  to  H2  He,,  continues  as  N 
increases  remains  to  be  seen.  In  the  bulk,  H2  penetration  is 
observed  but  the  extent  is  not  known.* 

Following  the  earlier  discussion  concerning  sampling 


efficiency  in  computing  pure  cluster  density  profiles,  one 
finds  that  the  efficiencies  for  computing  the  H,  and  He 
mixed  cluster^  density  profiles  scale  as  N  ~ 2  and 
(N  —  1  )N respectively.  This  fact  is  manifested  by  the 
relatively  large  fluctuations  in  the  H:  density  profile  as  N 
becomes  large,  cf.  Fig.  7.  For  all  mixed  clusters  small  bin 
sizes  were  employed.  A#  <0  03  A.  These  smaller  bin  sizes 
yield  a  more  detailed  picture  of  density  profiles  but  tend  to 
give  larger  statistical  error,  especially  near  the  origin.  How¬ 
ever,  despite  both  this  and  the  relative  inefficiency  of  binning 
for  a  single  H2  panicle  in  an  N-particle  cluster,  the  statistical 
error  in  the  H2  density  profiles  is  quite  reasonable  to  within 
about  1  A  of  the  cluster  center.  Even  for  the  largest  cluster, 
H2  He,, ,  computed  statistical  error  in  each  bin  is  10%  or  less 
from  l  to  10  A.  (For  He,  the  precision  is  of  course  much 
better:  less  than  5%  from  0.5  to  12  A  for  H2He„.)  Once 
again,  we  point  out  that  computational  cost  is  not  large,  e  g.. 
27  Cray/X-MP14  minutes  for  the  H,  He,,  calculation. 

Figures  6  and  7  demonstrate  that  H,  is  delocalized 
throughout  the  cluster,  as  is  the  case  for  He  in  both  the  pure 
and  mixed  clusters.  This  s  not  surprising  considering  the 
lighter  mass  of  H,  and  the  similarity  of  the  H,  -He  and  He- 
He  potentials.  Perhaps  most  interesting  is  that  the  greatest 
probability  is  not  observed  at  the  origin  but  peaks  at  some 
distance  from  the  center.  This  point  of  maximum  H2  proba¬ 
bility  density  increases  from  roughly  1  to  5  A  as  the  cluster 
increases  from  5  to  20  particles,  paralleling  the  increase  in 
cluster  size.  Note  that  for  H2  He2  the  probability  density  for 
H:  simply  mimics  that  of  He  (Fig.  6(a)  ]  and  that  both  the 
H,  and  the  He  density  profiles  are  similar  to  that  of  He, 

( F  |.  4) .  This  shows  that  H2  can  take  the  place  of  a  He  in  any 
position,  including  tne  collinear  arrangements  which  lead  to 
the  unique  density  profiles  we  have  observed. 

Table  V  summarizes  a  quantitative  analysis  of  the  H; 
probability  peak.  The  second  column  measures  the  absolute 
position  of  the  peak,  RPm  (relative  toR, )  and  clearly  reflects 
the  increase  in  size  of  the  cluster  as  N  increases.  The  third 
column  measures  the  degree  to  which  the  peak  in  probability 
manifests  itself:  pm/p0  =  p»,  (RPm)/pHl  (R  *  0)  is  the  ratio 
of  the  maximum  H2  probability  density  to  its  value  at  the 
origin.  ( For  H2  He,  the  second  maximum  is  used  as  /?„  ) 
Given  the  statistical  errors  in  p  near  the  origin  and  the 
widths  of  the  observed  peaks,  the  quantities  in  Table  V  are 
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estimated  with  some  degree  of  uncertainty,  i.e.,  5%-10%. 
However,  allowing  for  the  qualitative  nature  of  the  results,  it 
is  obvious  that  the  peak  at  Rtm  becomes  more  pronounced 
for  the  larger  clusters,  N>7.  Among  the  larger  clusters,  the 
Hj  He, ,  peak  value  of pm  /p0  is  significantly  greater  than  the 
others,  as  is  evident  also  from  Fig.  7.  The  fourth  column 
gives  the  ratio  of  the  He  density  at  RPm  to  its  value  at  the 
origin,  i.e.,pM.  (RPm  )/pH,  (R  =0).  The  fifth  column  shows 
the  amount  of  helium  contained  inside  the  peak  position,  i.e.. 


FIG.  6.  Normalized  He  and  H,  density  profiles  for  the  mixed  clusters  1 1 1 
H,  He, ,  ( b )  H,  He, .  ( c )  H,  He, .  Dashed  and  solid  curves  refer  to  He  and 
H,,  respectively,  and  are  normalized  to  one  particle  per  cluster  The  bin 
sizes  are  <0.03  A  and  R  is  the  distance  from  the  (cometnc  center  of  tbe 
cluster. 


S0'm4irpH'(R)R2dR.  Omitting  the  special  H3He,  case, 
these  last  two  columns  of  results  show  that  the  H:  molecule 
tends  to  reside  at  greater  distances  from  the  cluster  center  for 
the  larger  clusters.  However,  the  H,  molecule  shows  no  sig¬ 
nificant  tendency  to  move  further  towards  the  duster  ed  sc  as 
the  size  increases  from  13  to  19  He  atoms;  the  peak  becomes 
less  pronounced  (third  column.  Table  V)  and  its  position 
relative  to  the  helium  density  only  increases  by  a  email 
amount  (last  two  columns).  In  addition,  the  position  of  'hr 
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FIG.  7.  Normalized  He  and  H,  density  profiles  for  the  mixed  dusters:  ( a ) 
H,He,..  (b)  H,  He, , ,  (c>  H.He„  Dashed  and  solid  curves  refer  to  He  and 
Hj.  respectively.  For  further  details  see  Fig.  6. 


TABLE  IV  Size  comparisons  between  pure  and  mixed  clusters.’ 


Cluster 

(<**>) 1/1 

Cluster 

He 

H, 

(( *,))'n 

P 

HR1))'’1 

P 

He, 

6.1* 

Hi  He, 

5.67 

0417 

6.05 

0.332 

He. 

5.44 

Hi  He, 

5.41 

0.417 

5.81 

0.526 

He, 

4  98 

H,He. 

4.7J 

0.41S 

5.17 

0.570 

He. 

5.11 

Hi  He, 

5.02 

0.490 

5.37 

0.579 

He,, 

5  40 

H,H*i, 

5  39 

0.492 

3.77 

0.590 

He,, 

5*5 

H,He„ 

541 

0.494 

6.12 

0.657 

He*, 

5.93 

H,He„ 

596 

0.494 

6.34 

0601 

■  The  quantity  ((K1))1'1  is  the  rms  distance  from  the  clutter  center,  in  A.  P  is  the  cumulative  probability  of 
finding  the  particle  beyond  the  spherical  shell  in  which  half  the  part  idea  are  found  on  average  ( see  text  I . 
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TABLE  V.  H.  probability  peak  characteristics.* 


Cluster 

P~/Po 

(H,l 

P~/p« 

IHe! 

%  He 

inside 

H.He, 

2.3 

0.86 

0.76 

9  2 

H.He, 

1.3 

1.05 

0.89 

2.1 

H.He, 

27 

1  16 

0.38 

16.9 

H.He. 

2.9 

1  28 

0.83 

16.7 

H.He,, 

40 

I  32 

0.67 

29.5 

H.He,, 

4.4 

l  68 

0.65 

35  8 

H.He,, 

5  1 

1.38 

0.61 

40.3 

R„_  is  the  position  of  the  H.  maximum,  in  A.  pm/p„  is  the  ratio  of  the  peak 
probability  to  the  central  ( R  =  0)  probability,  for  H.  and  He.  respective¬ 
ly  The  %  He  inside  is  measured  with  respect  to 


peak  lies  well  inside  the  cluster  surface,  i.e.,  at  least  6 0%  of 
the  integrated  He  density  lies  beyond  the  H,  peak  There¬ 
fore,  the  Hj  density  profiles  show  some  localization,  most 
noticeable  for  the  larger  clusters  ( >  1 2  He  atoms ) ,  and  indi¬ 
cate  absorption  between  the  cluster  center  and  its  extenor 
This  perhaps  reflects  the  fact  that  H,  has  been  found  to  pene¬ 
trate  bulk  liquid  helium.' 

Figures  8  and  9  show  the  panicle  separation  probability 
density  functions  p(r)  for  the  He-He  (dashed  lines)  and 
H2-He  separations  (solid  lines).  These  show  the  same  gen¬ 
eral  trends  with  size  which  are  apparent  in  Figs.  4  and  5 
Differentiation  between  H2  and  He  is  not  large.  However, 
the  H2-He  distnbutions  are  consistently  smaller  at  smail 
separations  and  consistently  more  diffuse  at  large  separa- 


r  (A) 


FIG  8  Particle  separation  probability  density  function  pi ')  for  the  mn« 
clusters  la)  H.He,,  (b)  H.  He,,  (c)  H,He,  Solid  lines  refer  to  H  -He 
dished  lines  to  He-He.  The  normalization  and  bin  sizes  are  as  in  F:g  a 
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FIG  9  Particle  separation  probability  density  function  pi  r)  for  the  mutd 
clusters:  (al  H.  He,,.  (b)  H.He,,,  (cl  H.  He,.  Solid  lines  refer  toH.-He. 
dashed  lines  to  He-He.  The  normalization  and  bin  sizes  are  as  in  Fig  4 


5  10  15 

r  (A) 


tions  than  the  He-He  distributions.  The  trend  at  small  r  is  the  maximum  probability  of  finding  H:  is  somewhat  re- 
easily  understood  from  the  interaction  potentials,  cf.  Fig.  I,  moved  from  the  cluster  center. 

where  the  H,  -He  potential  becomes  highly  repulsive  sooner  As  is  the  case  for  the  density  profiles,  the  particle  separa- 

than  that  of  He-He.  The  behavior  at  large  r  is  less  transpar-  tion  distributions  p  stand  out  for  Hj  He,  j .  The  greater  rela- 

ent.  The  H2-He  potential  has  a  greater  well  depth  and  is  tive  diffuseness  of  p  for  H,-He  vs  He-He.  clearly  visible  in 

more  attractive  at  larger  separations.  On  the  other  hand,  the  Fig.  9,  is  most  pronounced  for  H :  He,,  amongst  the  larger 

lighter  mass  of  H2  tends  to  reduce  H2  -He  binding  vs  He-He  clusters.  The  implication  here,  as  with  the  density  profiles,  is 

binding  This  factor,  as  well  as  the  fact  that  helium  atoms  can  that  Heu  possesses  a  somewhat  enhanced  degree  of  stability 

get  closer  to  themselves  than  to  H,.  may  also  explain  why  with  respect  to  H,  penetration. 


S~\  H, He,, 

0  *  5 

_  H.He 

/  \  2  12 

0.12- 

jf'K  2  19 

I  \ 

0.09- 

_ °/<C 

(  \ 

J  \ 

C^-Q.06  - 

j  \ 

I  \ 

i - r 
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1  X 
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TABLE  VI.  Wave-function  parameters  for  He,  (m  atomic  units). 


N 

3 

4 

3 

7 

13 

14 

20 

a!  10 

-0.102  29 

-  0. 102  29 

-0.1*3  63 

-0.107  00 

-  1  059  5 

-  1  070  6 

-  0  600  74 

b 

-  1.3402* 

-  1.340  2* 

-  1  306  76 

-  1.4105* 

-  0  975  8* 

-  1  046  66 

-  1  046  64 

a 

1  183  70 

1.183  7$ 

1.035  09 

0  526  91 

0  603  14 

0  545  00 

0  545  03 

0.146  30 

0.146  30 

0  146  30 

-  0  971  90 

-  1  459  87 

-  1  323  61 

-  1  308  01 

<2, 

-  20.0*6 

-  23.447 

-  23.034 

-  23.037 

-  29  939  9 

-  30  193  5 

-  38  864  6 

J. 

131  69 

139  923 

158.499 

158.50 

200  994 

199  831 

310  061 

a, 

-  746  490 

-  742.004 

-  7*2.170 

-  782.170 

-  851  450 

-  845  506 

-  1  370  01 

<*4 

1  709.00 

1  710  60 

1  *42.10 

1  *42.10 

1  809  01 

1  808  9* 

2  484  45 

-4  404  | 

-4  403.6 

-  4  344.3 

-  4  344  3 

-  4  352.3 

-  4  354  1 

-  3  674  0 

IV.  CONCLUSIONS 

This  first  quantum-mechanical  study  of  an  impunty 
species  in  a  rare-gas  cluster  shows  the  lower-mass  Hj  impu¬ 
rity  to  be  extensively  delocalized  in  finite  helium  clusters. 
Considering  the  similar  mass  and  interaction  potential  with 
He  (cf.  Fig.  1 ),  this  is  consistent  with  observations  of  !He 
delocalization  in  bulk  helium  (4He).24  However,  much  less 
is  known  about  H2  or  any  other  molecular  impunty  in  bulk 
helium.  It  will  be  useful  to  follow  this  trend  for  cluster  sizes 
larger  than  those  studied  here,  to  ascertain  whether  the  H, 
species  is  surface  attached,  or  embedded  in  the  interior,  as 
appears  to  be  the  case  thus  far.  z'  we  approach  bulk  systems. 
In  addition,  possible  localization  v  "the  H2  and  of  other  im¬ 
purities  at  larger  cluster  sizes  is  of  interest.  Finally,  incorpor¬ 
ating  the  weak  anisotropy  of  the  Hj-He  potential  is  worth 
considenng. 

The  wave-function  forms  employed  here  for  both  the 
mixed  clusters  H,  -HeN  and  the  pure  HeN  reference  systems 
present  a  radical  departure  from  previous  analytic  forms 
used  in  variational  Monte  Carlo  ground-state  studies.  No 


one-particle  term  is  employed,  and  the  two-particle  term  is 
explicitly  constructed  by  fitting  to  the  bound  numerical  solu¬ 
tion  of  HjHe,  which  allows  much  greater  sensitivity  to  the 
long-range  part  of  the  potential  than  in  previous  forms.  Ac¬ 
curate  representation  of  the  diffuse  tail  of  the  two-body  fac¬ 
tors  appears  to  be  extremely  important  in  obtaining  accurate 
energies.  We  note  here  that  for  He,  and  He,  the  energies  he 
below  the  previously  computed  exact  GFMC  values.  In  Ref 
1 1,  we  find  that  this  results  from  the  use  by  GFMC  of  the 
1979  potential  which  is  slightly  shallower  than  the  most  re¬ 
cent  one  ( 1987)  employed  here.  The  high  level  of  accuracy 
obtained  here  at  the  VMC  level  with  just  the  two-particle 
term,  even  for  N  =  20,  also  suggests  that  addition  of  the 
three-particle  terms  known  to  be  important  for  the  larger 
clusters  may  well  provide  a  new  level  of  accuracy  at  the 
VMC  level  for  these  larger  systems  as  well.  It  will  be  very 
interesting  to  see  how  well  the  lack  of  one-particle  term  does 
at  larger  N.  This  term  has  been  required  in  previous  studies 
by  us  and  others  because  of  the  modified  liquid-state  asymp¬ 
totic  form  imposed  on  the  two-particle  term.  Note  that  the 
new  form  employed  here  does  have  the  flexibility  to  revert  to 


TABLE  VII.  Wave-function  parameters  for  H,  He,  (in  atomic  units). 


He-He  parameters  for  H,  He, 

N 

2 

3 

4 

6 

12 

13 

19 

a/ 10 

-  0.577  43 

-  0.083  449 

-0.391  86 

-  0.467  82 

-  1  059  5 

-  1  024  1 

-  0  840  34 

b 

-  1  078  43 

-  1  335  66 

-  1.590  36 

-  1.372  76 

-  1  016  85 

-  1023  28 

-0  976  41 

a 

0  970  574 

1.171  530 

0.791  906 

0  645  7*9 

0  603  138 

0  60*  093 

0  553  31  ? 

Jr, 

-  I  308  32 

0.114  79 

0.638  34 

-0.758  00 

-  1.508  01 

-  1  *05  82 

-  1  386  08 

a, 

-  23.932 

-  20.046 

-  25.642 

-  24  427 

-  30.753 

-  33  700 

-  38  875 

a: 

103.288 

151.690 

157,665 

158  838 

202.829 

226.526 

310  059 

j, 

-  734.485 

-  746.490 

-  782.390 

-  782  0*9 

-  852.536 

-  307.670 

-  1307  01 

jt 

1  715.65 

1  709.00 

1  842.20 

1  842.13 

1  813.57 

I  792.85 

2  4*4  45 

0, 

-4*01.3 

-4  404.1 

-4050.3 

-4  040.3 

-  4  348.4 

-  *  350.0 

-  3  674  6 

H,  -He  parameters  for  H,  He, 

N 

2 

3 

4 

6 

12 

13 

19 

<3/ 100 

-  0.657  03 

-0.050*60 

-  0.050  865 

-0.362  *3 

-0  362  *3 

-0411  31 

-  0  286  01 

b 

-  0.625  84 

-  1  433  77 

-  1.411  75 

-  1.120  54 

-0.730  13 

-0.715  35 

-  0  685  77 

a 

1.400  470 

1.147  770 

1. 133  880 

1.151  390 

1  151  390 

1.153  680 

1  069  830 

a. 

-  2.139  62 

-0  657  03 

-0.164  98 

-0  879  13 

-  0  927  28 

-0.125  0* 

-0  903  15 

<7, 

-6  163 

-  14.0*3 

-  15.359 

-  14  923 

-  12.744 

-  14  607 

-  14612 

112.344 

121.898 

121.843 

123.210 

122.727 

131.758 

131  757 

«i 

—  719. 100 

-  710.496 

-710.504 

-  710  000 

-  711.061 

-  730.737 

_  730  73- 

<7« 

1  723.29 

1  726  9* 

l  726.97 

1  727  10 

1  726.35 

1  715.99 

1  715  99 

a. 

-  7  455.6 

-  7  454.4 

-  7  454  2 

-  7  454  4 

-  7  454.7 

-  7  45*  0 

-  7  458  0 
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The  H,-4He  dimer  and  small  4He  clusters  are  studied  using  Monte  Carlo  sampling  techniques.  Wa 
consider  alternative  wave-function  forms  in  order  to  obtain  high  accuracy  efficiently.  For  the  «m»Hg 
systems-  both  guided  and  unguided  Metropolis  walks  are  used  and  the  efficiencies  are  studied.  Of  partic¬ 
ular  concern  is  accurate  sampling  at  small  particle  separations  and  the  behavior  of  the  local  energy  in 
this  regime.  As  a  final  step,  we  compute  exact  energies  by  a  diffusion  Monte  Carlo  method.  We  obtain 
converged  energies  significantly  below  the  Creen's-function  Monte  Carlo  values,  which  employed  an  ear¬ 
lier  He-He  potential  with  a  slightly  shallower  well.  For  He>  and  Hejo,  the  Green’s-function  Monte  Carlo 
energies  are  reproduced  when  employing  the  same  potential.  However,  for  the  112-atom  cluster,  our 
converged  energy  lies  below  the  Green's-function  Monte  Carto  value.  Second-order  estimates  of  the  ex¬ 
act  density  profiles  and  particle  separation  distributions,  p,  are  also  determined.  For  the  14-  and  20-atom 
clusters,  second-order  estimates  of  p  show  enhanced  structure  in  comparison  to  variational  Monte  Carlo 
results.  Statistically  meaningful  oscillationa  in  the  second-order  estimates  of  the  exact  density  profiles 
are  not  observed. 

PACS  numberts):  36.40.  +d.  67.40.Db.  02.70. -c,  03.65.G* 


,  I,  INTRODUCTION 

Atomic  and  molecular  clusters  have  become  of  interest 
to  both  theorists  and  experimentalists  [1],  Of  particular 
concern  are  structure,  phase  transitions,  binding  energies, 
and  excitation  spectra,  and  the  behavior  of  these  proper¬ 
ties  as  the  bulk  is  approached. 

We  are  interested  in  studying  atomic  and  molecular 
clusters,  both  pure  and  with  impurities  attached,  using 
Monte  Carlo  techniques.  Such  approaches  thus  far  pos¬ 
sess  the  greatest  possibility  of  yielding  high  accuracy  for 
theoretical  methods.  To  enhance  the  capabilities  of 
Monte  Carlo  for  these  systems,  we  consider  alternative 
wave-function  forms  and  the  efficient  optimization  of 
wave-function  parameters  in  studying  weakly  bound 
quantum  clusters.  To  start  with,  we  study  the  HtHe 
complex  (from  here  on  4He  is  assumed).  This  system  is 
quite  useful  as  it  provides  a  very  weakly  bound,  highly 
repulsive  potential,  test  case  for  the  initial  wave-function 
form  we  employ. 

As  a  further  development,  we  employ  diffusion  Monte 
Carlo  (DMC).  We  use  this  approach  to  compute  exact 
ground-state  energies  for  helium  clusters  with  the  most 
up-to-date  potential.  In  addition  to  increased  accuracy  in 
the  energy  and  structural  features  such  as  the  density 
profile,  the  DMC  approach  serves  to  provide  benchmarks 
for  evaluating  wave-function  quality.  This  is  pertinent 
for  the  helium  clusters  for  which  exact  energies  resulting 
from  the  most  recent  pair  potential  have  not  yet  been 
computed. 

The  remainder  of  the  paper  is  organized  as  follows.  In. 
Sec.  II  we  discuss  Monte  Carlo  integration  techniques, 
and  in  Sec.  Ill  the  exact  diffusion  Monte  Carlo  approach. 
The  wave-function  forms  and  optimization  technique  we 
employ  are  presented  in  Sec.  IV.  In  Sec.  V  we  present  re¬ 
sults  for  a  range  of  small  clusters  ( /V  120  and  N  =  112). 
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Section  VI  presents  conclusions  concerning  the  wave- 
function  accuracy  and  sampling  efficiency  for  a.V  tech¬ 
niques. 

II.  MONTE  CARLO 
INTEGRATION  TECHNIQUES 

Multidimensional  integration  is  performed  by  Monte 
Carlo  in  order  to  obtain  wave-function  expectation 
values.  This  is  achieved  by  sampling  points, 
R=fr,,rj, . . .  ,rM),  from  a  probability  density  function 
p(R).  Expectation  values  of  coordinate  operators  .4  (R) 
are  then  computed  as 

_  M 

AM=ATl  2  AM,) ,  (1) 

i“l 

with  { R,  j  sampled  from  p.  As  M  becomes  large,  AM  ap¬ 
proaches  the  average  of  A  (R)  over  p(R). 

We  employ  two  variants  of  the  Metropolis  walk  [2,3] 
to  sample  p.  The  first  of  these  is  the  widely  used  and  very 
simple  “unguided”  walk.  For  a  point  at  R,  a  new  point  is 
sampled  from  a  transition  probability  density  r(R— »R’) 
which  is  simply  constant  within  a  cube  and  zero  outside. 
Thus  moves  are  g:  en  by 

R'=R  +  A(f— 0.5!) ,  (2) 

where  $  is  a  vector  whose  components  are  uniform  ran¬ 
dom  variates  between  0  and  I. 

This  "unguided”  walk  attempts  to  move  uniformly 
through  coordinate  space  without  regard  to  the  form  of 
|'F|2.  Therefore  a  more  efficient  scheme  of  choosing  at¬ 
tempted  moves  is  likely.  This  is  the  basis  of  a  guided  or 
“smart'’  Metropolis  walk,  which  is  also  known  as  impor¬ 
tance  sampling.  We  now  choose  the  transition  density  to 
be 
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r{R-.R')=(4«-A)",#/,expj  -[R'-R-AF(R)]V4A|  . 

(3) 

Sampling  from  this  transition  density  requires  that 

R'=R+AF(R)+v''2A;r .  (4) 


probability  density  function  /  ==*<*»,  where  *  is  a  tnai 
wave  function.  Rewriting  Eq.  (5)  in  terms  of  /  gives 

-|f  =  -^^/+DF-[/F(R)]  +  {£t(R)-£J,l/,  (7) 

where 


The  components  of  X  are  Gaussian  random  numbers  with 
a  mean  of  zero  and  a  variance  of  unity,  which  we  obtain 
by  the  Box-Mueller  algorithm  [4].  The  “guiding  force,” 
F= V|'p|2,  acts  to  push  moves  in  the  direction  of  most 
rapidly  increasing  |4*|2. 

The  major  consideration  for  the  approaches  discussed 
here  is  the  value  of  A  which  yields  the  most  efficient  sam¬ 
pling.  The  optimum  choice  lies  between  a  small  value, 
which  yields  a  high  acceptance  rate  but  a  large  degree  of 
correlation  between  moves,  and  a  large  value,  which  gives 
large  attempted  displacements  but  a  very  small  accep¬ 
tance  rate  so  that  correlation  is  large.  The  optimum  A  is 
often  taken  to  be  that  which  yields  an  average  acceptance 
rate  (or  ratio)  of  roughly  0.5.  Here  we  consider  a  further 
quantitative  measurement  of  the  efficiency  with  which 
configuration  space  is  sampled,  namely,  the  average  dis¬ 
placement  of  moves  during  the  walk,  (A  R  ).  (A  rejected 
move  contributes  a  value  of  zero  to  this  average.) 


III.  DIFFUSION  MONTE  CARLO 

Although  the  Metropolis  algorithm  provides  a  means 
for  computing  expectation  values  of  a  given  wave  func¬ 
tion,  accuracy  is  limited  by  the  quality  of  However, 
exact  Monte  Carlo  approaches  are  well  known.  These 
approaches  are  often  genetically  referred  to  as  quantum 
Monte  Carlo  and  fall  into  two  categories.  Green’s* 
function  Monte  Carlo  (GFMC)  and  diffusion  Monte  Car¬ 
lo.  The  former  has  been  applied  to  a  wide  range  of  prob¬ 
lems  and  derives  from  consideration  of  the  time- 
independent  Schrodinger  equation.  Initial  work  was  on 
the  helium  atom  [5]  and  liquid  helium  [6,7],  and  later  ap¬ 
plications  include  electronic  structure  calculations  [8,9] 
and  Computations  on  helium  clusters  [10-12]. 

DMC  starts  with  the  time-dependent  Schrodinger 
equation  in  imaginary  time  and  has  been  employed  most¬ 
ly  in  electronic  structure  calculations  [13-17].  Recent 
work  has  also  included  helium  clusters  and  other  van  der 
Waals  species  [18-20].  The  DMC  approach  we  employ 
is  very  similar  to  that  of  Reynolds  el  al.  [15]  and  is  out¬ 
lined  below. 

Writing  the  Schrodinger  equation  in  imaginary  time, 
i — t  //,  and  setting  =  1  we  have 


8d»(R ,/) 
3f 


=<//-£*  M>(R,f) . 


(5) 


The  reference  energy  EK  only  affects  the  (imaginary)  time 
dependence  of  <t>(R,r).  It  is  easily  shown  that  at  large  t 
the  ground  state  dominates,  leaving 


<t>=exp[-«E0 -£„>]*„. 


(6) 


Note  that  the  choice  of  £*  =  E0  is  useful  in  removing  the 
time  dependence  from  the  asymptotic  solution. 
Importance  sampling  is  implemented  by  choosing  a 


//=-DVJ+F,  D72hX <2mir,V21 

i 

El='¥~'H\ V  is  the  local  energy,  and  once  again 
F= Vln|'F|2.  Note  that  terms  on  the  right-hand  side  of 
Eq.  (7)  correspond  to  diffusion,  drift,  and  branching,  re¬ 
spectively.  The  asymptotic  form  of  /  follows  from  Eq.  (6) 
and  it 

/(R,r)=exp[-f(£’0-EJ,)]'P(R)A0(R)  .  (8) 

When  /  takes  this  form,  expectation  values  over  /  are  in¬ 
dependent  of  t,  i.e., 

(A)f  =  J* f(R,t)A  (R)dR  J //(R,rWR 

=  / WRJdolRMfRldR/  f  ♦(R)d0(Rk/R  . 

For  /t  (R)  =  £t(R),  it  is  easily  shown  that  ( Et  )f~EQ, 
so  that  the  ground-state  energy  is  obtained  as  the  average 
of  El  over  /.  The  time  development  of  /  is  given  by 

/(R\r+r)=  J dR/(R,r)G(R-*R\r) ,  (9) 

where  the  Green's  function  G  describes  a  move  from  R  to 
R’  in  time  r.  The  Green’s  function  is  a  solution  of  Eq.  (7) 
with  the  boundary  condition  G(R— R',0)=6(R-R*). 
For  all  but  a  few  simple  Hamiltonians,  the  Green's  func¬ 
tion  is  unknown.  Here,  we  employ  an  analytic  ‘’short- 
time”  approximation  to  G  which  takes  the  form 

Ga(R— R*,r)=(4ffDT)-J#/2 


Xexpj  —  [R‘  — R  — DrF(R)]2/4Dr] 
Xexp(-r[[£t(R’)  +  £t(R)]/2-£<})  . 


(10) 

The  approximate  nature  of  Ga  is  clear  from  Eq.  (10):  dur¬ 
ing  the  course  of  a  move  from  R  to  R'  in  time  r,  the  drift 
[determined  by  F(R)]  and  branching  (dependent  on  £t) 
are  assumed  to  be  constant.  While  error  in  G,  vanishes 
as  r— 0  [21-23],  for  the  practical  case,  i.e.,  r=?*=0,  the 
asymptotic  /  only  approximates  ♦do,  and  computed  re¬ 
sults  will  differ  from  the  corresponding  averages  over 
'PSq.  This  difference,  referred  to  as  time-step  bias,  may 
either  be  removed  by  extrapolation  or  made  insignificant 
by  using  values  of  r  such  that  the  bias  is  less  than  the  sta¬ 
tistical  error. 

To  reduce  time-step  bias,  an  acceptance-or-rejection 
step  is  employed  [15].  That  is,  moves  are  accepted  with  a 
probability  A  given  by 

/4(R  —  R',r)  =  minj  l,u>(R  —  R',r)|  ,  ill) 


and 


u’(R  — R',t)  = 


jW  R');|GU(R*  — R,r) 
|«MR)|-’G,(R  —  R'.t) 


i:i 
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Including  the  factor  A  plays  an  important  role.  As  'l'  ap¬ 
proaches  the  exact  solution  ^0,  the  branching  becomes 
constant  and  Ga  is  essentially  the  transition  density  given 
in  Eq.  (3)  with  A  =  /)r.  In  this  instance,  DMC  reduces  to 
a  guided  Metropolis  walk  and  f(=\<t>Q\1)  is  sampled 
without  time-step  bias — for  any  value  of  t.  Use  of  A  has 
been  found  to  greatly  reduce  time-step  bias  [24]  because 
the  acceptance-or-rejection  step  eliminates  time-step  bias 
to  the  extent  that  ♦  resembles  ^q. 

We  conclude  this  section  with  a  discussion  of  several 
technical  details.  In  the  DMC  [and  Metropolis-walk  or 
variational  Monte  Carlo  (VMC)]  computations,  ten  in¬ 
dependent  ensembles  of  100  walkers  are  propagated  in 
parallel  yielding  ten  independent  estimates,  of,  for  exam¬ 
ple,  the  energy,  from  which  the  average  and  its  statistical 
error  are  obtained.  If  greater  precision  is  desired,  more 
runs  are  performed  in  this  manner.  This  structure  is  use¬ 
ful  in  that  for  each  set  of  runs  an  estimate  of  statistical  er¬ 
ror  unbiased  by  serial  correlation  is  computed. 

Updating  the  reference  energy  ER  is  useful  for  minim¬ 
izing  the  time  dependence  of  the  ensemble.  For  a  given 
number  [or  population,  /MO)]  of  points  sampled  from 
/ (0)  =  at  the  beginning  of  the  DMC  walk,  we  have 

P(0)—  f  f(R,0)dR  .  (13) 


ing.  As  seen  from  Eq.  (10),  the  branching  factor  for  a 
move  from  R  to  R'  is 

*(R.R')  =  exp(-r|[£t(R')  +  £t(R)]/2-£J,})  .  (16) 

Branching  may  be  implemented  by  obtaining  an  integer 
.V/  =  intji>(R,R')  +  £],  where  int(x)  is  the  largest  integer 
that  is  <  x  and  where  £  is  a  uniform  random  variate  uni¬ 
formly  distributed  between  0  and  1  so  that  M  ~b  on 
average,  and  creating  M  copies  of  walkers  at  R'.  Alterna¬ 
tively,  one  may  assign  a  weight  a>(R')  =  b<R,R')iu(R)  to 
the  walker  at  R’.  Since  Si  equals  6  only  on  average,  as¬ 
signing  (or  carrying)  wetghts  would  seem  preferable. 
However,  earned  weights  diverge  towards  0  or  <»=  as  the 
walk  proceeds,  giving  rise  to  the  possibility  that  an  en¬ 
semble  may  he  effectively  composed  of  a  few  walkers  with 
very  large  relative  weights.  In  this  event,  the  sampling  is 
inefficient  as  only  a  few  of  the  many  walkers  being  pro¬ 
pagated  contnbute  to  the  computed  averages.  Therefore 
we  employ  a  combination  of  carrying  weights  and  copy¬ 
ing.  Weights  are  earned  unless  w  £  u;min  or  w>  If 
either  of  the  bounds  are  exceeded  then  Sfw  —  mt(  w  +  £ ) 
copies  are  made,  and  for  w  >  u’m„  each  copy  is  assigned  a 
weight  of  w/Mw.  For  the  DMC  results  reported  here, 
“’nut  ~0-  *  -0.4  and  — 2. 


From  the  asymptotic  form  of  /  it  is  easily  seen  that 

f>(r)=exp[-r(£0-£'J,)]P<0)  .  (14) 

Note  that  Eq.  (14)  indicates  that  an  estimate  of  £0  may 
be  obtained  from  the  change  in  the  ensemble  size  over 
time.  This  estimate,  usually  referred  to  as  the  growth  en¬ 
ergy  (£c),  often  possesses  a  different  dependence  on  the 
time  step  than  does  the  average  of  the  local  energy  EL . 
Therefore,  to  reduce  the  long-term  growth  or  decay  of 
ensemble  size,  at  each  time  step  we  perform  a  short  run 
to  estimate  £c  and  then  set  ER  equal  to  this  estimate 
when  computing  the  reported  results. 

Another  point  concerns  the  renormalization  of  the  en¬ 
semble  population  P.  Even  when  ER  is  equal  to  the 
growth  estimate  of  £0,  fluctuations  in  P  arise  from  fluc¬ 
tuations  in  Ec,  which  are  in  turn  caused  by  variations  in 
the  local  energy,  as  the  ensemble  is  propagated.  If  the 
statistical  error  in  £G  over  the  ensemble  is  <7^,  then 
from  Eq.  (14)  the  relative  statistical  error  in  P  is  seen  to 
increase  proportionally  to  time  as 

<7 p/P  —  ta E{,  .  (15) 

Therefore,  in  keeping  the  ensemble  size  reasonable,  it  is 
useful  to  renormalize  the  population  [back  to  P0  (  =  100)] 
at  intervals  during  the  walk.  However,  as  noted  previ¬ 
ously  [25],  renormalization  introduces  an  error  which  de¬ 
creases  as  the  frequency  of  renormalization  decreases. 
(Generally,  this  error  is  not  noticeable  unless  the  propa¬ 
gation  time  between  renormalizations  is  very  short.) 
Here,  we  divide  each  run  into  blocks,  and  at  the  end  of 
each  block  the  population  is  renormalized  to  100  walkers. 
Wc  vary  block  propagation  times  ( 1 0*  —  1 0*  hartree"’)  to 
verify  that  “renormalization"  bias  is  negligible. 

The  final  point  concerns  the  implementation  of  branch¬ 


IV.  TRIAL  FUNCTION  FORM 
AND  OPTIMIZATION 


A.  Trial  functions 


We  seek  ground-state  wave  functions  for  bosonic  sys¬ 
tems.  Such  wave  functions  are  nodeless  and  therefore 
may  be  taken  as  everywhere  positive.  A  convenient  rep¬ 
resentation  takes  the  form 


'KRl^exp 


2 


r,(Ri 


where  in  the  completely  general  case 


(17) 


T,( R)=  £  t,(r„rj,Tk, . . .)  .  (18) 

I,;,*,  •  •  • 

If  </<*<•••  > 

In  practice  /  is  taken  to  be  <  3  so  that  the  wave  function 
is  reasonably  compact.  Since  the  potential  is  given  by  the 
sum  of  pairwise  interactions,  we  omit  Tx  and  instead 
start  with  two-body  terms.  (One-body  terms  have  been 
employed  in  previous  studies  of  clusters  [10,26,27],  but 
are  not  necessary.)  We  then  add  on  three-body  terms  for 
increased  accuracy  if  desired.  This  reflects  the  fact  ihat 
two-body  effects  should  be  most  important,  especially  for 
weakly  bound  clusters,  followed  by  three-body  effects, 
etc. 

As  is  common  practice,  we  employ  a  two-body  term 
which  is  a  function  only  of  particle  separations,  i.e.. 


T2iR)~  2  2  MM/)-  <19) 

‘■j  >./ 

U  <  j)  U  <  p 

TTi is  term  is  both  translationally  and  rotntionally  invari¬ 
ant,  i.e..  Pc  ,,,  T2=0  and  Lr2=0,  as  required  for  (he 
ground-state  wave  function.  Given  the  importance  of 
Iwo-body  effects,  it  is  useful  to  explore  optimal  forms  for 
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t2 ■  In  this  vein,  H2He,  the  first  “cluster”  in  the  H2Hev 
series,  provides  an  excellent  test  case  as  a  weakly  bound 
species  with  an  interaction  potential  very  similar  to  that 
of  He-He.  Therefore  a  form  of  r2  which  accurately  de¬ 
scribes  H2He  may  prove  advantageous  for  the  helium 
clusters  as  well. 

The  initial  form  of  r2  is  motivated  by  our  studies  of 
H2He.  For  the  interaction  potential,  we  use  a  Lennard- 
Jones  plus  van  der  Waals  fit  (28]  to  the  ab  initio  data  of 
K0  computed  by  Meyer,  Hariharan,  and  Kutzelnigg  [29]. 
Since  our  potential  is  a  function  only  of  the  He  distance 
from  the  H,  center  of  mass,  the  Schrodinger  equation  for 
this  system  can  now  be  reduced  to  a  one-dimensional 
problem.  A  numerical  solution  for  the  ground  state,  do*1, 
is  obtained  with  the  finite  difference  algorithm  of  Schatz 
[30].  The  next  step  consists  of  considering  forms  for  'F 
and  fitting  them  to  do'1.  The  accuracy  of  both  the  wave- 
function  fit  and  the  computed  energy  expectation  value 
allow  an  assessment  of  the  quality  of  'F. 

In  determining  a  useful  analytic  form  for  *4',  we  treat 
the  long-  and  short-range  behaviors  separately.  That  is, 
takes  the  form 

Wr) ,  (20) 

where  s  and  /  denote  the  short-  and  long-range  functions. 


respectively.  The  short-range  form  is  chosen  to  be  con¬ 
stant  at  large  r.  Given  the  form  of  V  at  small  r,  a  natural 
choice  for  is 

4'J(r)  =  exp[/>(w)],  u  =r-1  , 

,  (21) 

P(u)  =  2  akuk  . 

k  -0 

The  bound  on  the  powers  included  in  P  results  from  the 
desire  to  limit  singularities  in  the  kinetic  energy  to  be  no 
greater  than  r~IJ,  given  that  this  is  the  dominant  singu¬ 
larity  in  our  potential  [28],  A  high-quality  form  of 
was  found  to  be 

V,  =  rbexp[ora]  .  122) 

Since  is  very  nearly  constant  at  large  r,  we  first  fit 
to  Infdo'1]  *n  this  domain  to  determine  a,  b,  and  a.  We 
then  determine  the  short-range  parameters  jak]  by 
fitting  in  the  highly  repulsive  and  in  the  well  regions. 
The  range  of  points  included  in  the  short-  and  long-range 
fits  determines  the  parameters,  which  are  listed  in  Table 
I.  We  find  that  the  analytic  wave  function  ‘F  is  nearly  in¬ 
discernible  from  the  numerical  solution  do'1,  and  the  en¬ 
ergy  is  reproduced  to  high  accuracy,  —0.02443  versus 
—  0.024  61  K — an  error  of  only  0.7%. 


TABLE  I.  Wave-function  parameters. 


Cluster 

Paramet«\^ 

H,He 

He,„ 

He20 

Henj(Tj) 

Hcm(7-j  +  r,> 

a 

-0.007  315  22 

-0.107061 

-0.062 

-0.01000 

-0.01400 

b 

-1.438143 

-1.04666 

-1.055 

-0.85000 

-0.8550 

a 

1.138  39 

0.559995 

0.545 

0.545031 

0.545031 

0Q 

0.1536 

-1.323  61 

-1.308  01 

— 1.30801 

-1.30801 

a  t 

-14.04255 

-30.1935 

-38.8646 

-38.8646 

-38.8646 

o 1 

121.49628 

199. 831 

310.061 

310.061 

310.061 

"3 

-710.897  88 

-  845.506 

-1370.01 

—  1 370.01 

-1370.01 

1726.975  6 

1808.98 

2484.45 

2484.45 

2484.45 

os 

-7454.4214 

-4354.11 

-3674.60 

-3674.60 

-3674.60 

Ao 

0.012 

0.012 

0.009 

a»o 

1.60 

1.60 

1.8 

'o 

5.0 

5.0 

4.5 

A, 

0.043 

0.043 

0.031 

<u, 

2.025 

2.025 

2.225 

3.225 

3.225 

3.225 

Cluster 


Parameters 

^  He, 

He, 

He;o 

aa 

-3.2319500 

-3.247  39 

-3.413  98 

a  i 

-0.068458  3 

-0.046985  0 

-0.019  1689 

a. 

8.55 

8.65 

9.20 

a. 

0.8400750 

0.914958 

0.873  622 

0. 

0.086  148  1 

0.069  567  6 

0.069  1040 

-632.979 

-682.409 

-675.000 

/, 

588.918 

590.085 

590.085 

ti 

-205.068 

-203.773 

-203.775 

h 

34.1321 

33.1709 

33.1633 

t < 

-2.628  22 

-2.46631 

-2.4S5  77 

r, 

0.060  1202 

0.053  3509 

0  058  8524 
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The  H:Hc  wave  function  gives  a  two-body  function, 

/2( r)  =  lnr  +  ar°+  2  akuk  ,  (23) 

k  -0 

which  we  have  employed  in  VMC  computations  of  heli¬ 
um  clusters  with  and  without  a  H2  impurity  [28].  The 
form  of  /2  is  structurally  similar  to  forms  used  previously 
[10,26,27],  The  differences  are  that  in  the  short-range 
form  all  powers  of  r-1  up  to  five  are  included  here,  and 
there  is  an  added  flexibility  in  the  long-range  form,  intro¬ 
duced  by  the  exponent  a.  As  discussed  above,  the  short- 
range  component  of  r2  is  based  on  our  H2-He  potential 
which  takes  a  Lennard-Jones  form  at  small  separations 
[28].  While  the  shape  of  this  potential  is  similar  to  that 
of  He-He  [31],  the  analytic  forms  are  quite  different. 
Therefore  it  is  of  interest  to  consider  entirely  new  forms 
for  r,  based  on  the  short-range  behavior  of  the  He-He  po¬ 
tential.  We  employ  here  the  most  recent  HFD-B(HE)  po¬ 
tential  of  Aziz,  McCourt,  and  Wong  [31]  for  all  calcula¬ 
tions  unless  it  is  explicitly  stated  that  the  earlier 
HFDHE2  [32]  potential  is  used.  . 

We  give  special  emphasis  to  regions  of  small  separation 
because  the»local  energy  generally  possesses  its  greatest 
fluctuations  as  r  becomes  small.  Therefore  a  two-body 
form  which  is  more  physically  correct  in  this  domain  will 
reduce  fluctuations  in  EL,  yielding  greater  precision,  and 
hopefully  greater  accuracy,  in  computed  energies.  Our 
new  r2  takes  the  form 

r2(r)=r(r)i»(r)  +  a0+a|r  +  y  lnr  ,  (24) 

where 

u(r)  =  atexp(  —  a„r—  /3„r2)  (25) 

mimics  the  short-range  form  of  the  He-He  potential  [31], 
and 

Tlr)=  2  tkrk  ,  (26) 

k  =o 

with  nk  =5.  We  have  chosen  y  =0  leaving  only  a  singu¬ 
larity  of  r  ~ 1  in  El,  arising  from  the  kinetic  energy.  Oth¬ 
er  permissible  values  are  y  <  }  which  give  an  r  ~2  singu¬ 
larity  but  a  finite  statistical  error  in  £t.  With  y  — 0,  r0a„ 
negative  and  large  gives  a  wave  function  which  is  very 
small  but  remains  nonzero  at  r  =0.  This  reflects  the  fact 
that  the  potential  converges  to  a  very  large  but  finite 
value  at  r  =0. 

Overall,  in  comparison  to  Eq.  (23),  f2  gives  added  em¬ 
phasis  to  domains  containing  small  particle  separations 
and  somewhat  less  emphasis  to  large  r.  It  is  hoped  that, 
by  directly  including  a  "potential-like"  function  in  'P,  the 
highly  repulsive  potential  term  will  be  better  canceled  by 
the  kinetic  contribution,  producing  smaller  fluctuations 
in  El  at  small  values  of  r. 

The  use  of  a  three-body  term  in  ground-state  wave 
functions  has  yielded  significant  improvements  in 
descriptions  of  both  the  liquid  [33-37]  and  clusters  [10J. 
Here,  we  employ  the  description  of  three-body  correla¬ 
tions  used  previously  in  microscopic  studies  of  quantum 
clusters  [10,26,27,38],  namely, 


r3(R)=Xo2  rjii)-  2 

'  !  t»»> 

2  nr.ml2-  2  |fo,)f,} 

»  j  (*i i 


with 

r ,</)=  2  SM,yv.  /=  0,1 . 

In  Eqs.  (27)  and  (28),  we  have 


|0(r)  =  (r  —  r0  )exp 


£,(r)  =  exp 


(27) 


(28) 


(29) 


Derivatives  of  T 3  are  evaluated  analytically.  This  is  actu¬ 
ally  faster  than  finite  difference  because  derivatives  by 
finite  difference  require  three  evaluations  of  the  exponen¬ 
tials  in  Eqs.  (29)  while  analytic  computation  requires  only 
one.  Overall,  adding  T j  to  F2  only  increased  computa¬ 
tion  time  by  about  85%  for  the  14-,  20-,  and  112-atom 
clusters. 


B.  Optimization 

Wave-function  parameters  are  optimized  by  hand  and 
by  conjugate-gradient  line  minimization.  Although 
crude,  varying  parameters  by  hand  is  quite  useful  in  com¬ 
plementing  more  sophisticated  techniques.  Since  initial 
wave-function  parameters  are  often  quite  poor,  hand  op¬ 
timization  can  quickly  yield  large  improvements.  The  re¬ 
sulting  wave  function  can  then  be  used  as  input  for  the 
more  sophisticated  optimization  techniques. 

The  conjugate-gradient  technique  [39]  seeks  a  local 
minimum  by  moving  in  directions  in  parameter  space 
which  are  conjugate  to  each  other,  leading  to  efficient 
convergence.  Essentially,  the  algorithm  consists  of  suc¬ 
cessive  line  minimizations  in  parameter  space.  This  pro¬ 
cedure  requires  the  computation  of  the  quantity  being 
minimized  and  of  its  first  derivatives  with  respect  to  the 
wave-function  parameters.  Here,  we  consider 


<*![£,■ -F,]2!'!') 

W) 


(30) 


This  quantity  is  useful  in  that  one  may  seek  either  a 
minimum  in  (he  energy  by  choosing  E%  «(EL)  or  in  the 
variance  by  choosing  £s  =  (£t>.  In  the  first  case, 
minimizing  the  energy  yields  global  accuracy,  and,  in  ibe 
second,  minimizing  fluctuations  in  the  local  energy  em¬ 
phasizes  local  accuracy  in  4'.  Since  computed  energies 
are  most  often  compared  in  discussions  of  accuracy,  we 
focus  on  minimizing  the  energy. 

The  integrations  in  Eq.  (30)  are  performed  by  averag¬ 
ing  over  a  fixed  vet  of  points  sampled  from  a  distribution 
'kyi2  corresponding  to  an  initial  set  of  parameters  [40]. 
Therefore  we  minimize  the  estimate  of  s2  given  by 
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2  [£t(Rl)-£|]2|4'(RJ/|4'0(Ri>lI 

*£  =  — - - -  .  (31) 

2  |4'(R,)/l'I'0(R()i2 

i-i 

The  major  consideration  for  the  stability  and  accuracy 
of  the  optimization  is  that  jy  accurately  approximate  s2. 
For  this  reason  the  ratio  i4'/*l'0|2  is  included,  reflecting 
the  fact  that  changes  as  parameters  are  varied,  al¬ 
though  this  requires  computing  parameter  derivatives  of 
>4>.  To  enhance  numerical  precision,  we  adjust  the  nor¬ 
malization  of  4'  so  that  2 i4V4'0|2  =  A/.  This  is  useful 
for  clusters  with  more  than  five  atoms  where  changes  in 
wave-function  parameters  have  a  large  effect  on  the  nor¬ 
malization  of  'V,  because  of  its  product  form,  and  could 
make  i4'/4/;)i2  uniformly  exceedingly  large  or  small.  The 
remaining  determinant  of  the  accuracy  of  sh  is  the  num¬ 
ber  of  points  in  the  fixed  sample.  While  a  large  number 
of  points  is  desired  for  high  accuracy,  M  is  limited  by 
considerations  of  computational  cost  and  memory  re¬ 
quirements.  M  is  chosen  .so  that  the  statistical  error  in 
the  average  of  EL  over  the  points  is  significantly  less  than 
the  desired  improvement  in  the  energy.  We  have  em¬ 
ployed  1000-2000  points  for  the  smaller  clusters,  He3_5  7, 
and  5000  points  for  He14  and  He,0. 

The  final  step  in  obtaining  reliable  optimizations  in¬ 
volves  updating  the  fixed  sample.  As  the  wave  function 
changes,  the  points  sampled  from  i4'0i2  become  a  poorer 
choice  for  computing  expectation  values  with  respect  to 
!4'l2.  This  is  most  noticeable  when  'Fq  is  of  poor  quality 
and  'F  changes  significantly  during  optimization.  One 
manifestation  of  this  degradation  of  the  fixed  sample  is 
divergence  of  the  energy  to  unrealistically  low  values. 
Therefore  we  have  found  it  useful  to  update  the  sample 
by  using  a  Metropolis  walk  to  generate  a  new  set  of 
points  sampled  from  the  current  wave  function  (which 
then  plays  the  role  of  ^0.)  Updating  is  implemented  after 
the  energy  has  converged  or  when  it  begins  to  diverge. 

While  the  conjugate-gradient  technique  has  been  suc¬ 
cessful  for  the  3-20-atom  clusters,  it  appears  to  be  much 
less  practical  for  larger  clusters.  (For  HeII2  we  started 
with  the  optimized  He20  parameters  and  only  reoptimized 
the  long-range  parameters  by  hand.)  As  the  number  of 
atoms  in  the  cluster  increases,  the  dimensionality  of 
configuration  space  which  must  be  represented  by  the 
fixed  sample  of  points  increases.  Therefore  larger  sam¬ 
ples  are  generally  required  for  the  larger  clusters, 

V.  RESULTS  AND  DISCUSSION 
A.  Metropolis  walks 

As  discussed  in  Sec.  II,  wave-function  expectation 
values  may  be  computed  by  the  Metropolis  approach. 
We  first  consider  H2He.  Since  the  wave  function  depends 
only  on  the  distance  of  He  from  the  midpoint  of  H2,  r, 
the  walk  takes  place  in  only  three  dimensions,  r. 

We  start  with  the  unguided  walk  employing  several 
values  of  A  in  the  range  11-40  A  and  M  =  5X107.  We 
find  that  sampling  l^l2  for  such  a  diffuse  system  is  not 


trivial.  The  energy  is  reproduced  reasonably  well,  the 
computed  value  is  generally  within  one  or  two  standard 
deviations  of  the  analytic,  and  the  average  error  is  0.1%. 
However,  we  find  sizable  errors  in  (r>  and 
rmt  =<  (r2)  )1/2-  The  smallest  value  of  A  produces  errors 
of  6%  and  7%  in  (r  >  and  rmi,  respectively.  Apparently, 
poor  efficiency  in  sampling  |4'|2  occurs  in  this  case. 
Upon  increasing  A  to  27  A,  giving  an  acceptance  ratio 
close  to  the  often  assumed  optimum  of  0.5,  errors  in  (r ) 
and  rrm  are  reduced  to  1%  and  2%,  respectively.  In¬ 
terestingly,  we  find  that  only  at  a  large  value  of  A,  40  A, 
and  a  relatively  small  acceptance  ratio  of  0.35,  is  accura¬ 
cy  in  <r )  and  rrm,  comparable  to  that  of  the  energy,  i.e., 
=0.1%. 

The  average  displacement  of  the  moves,  ( 1R  ),  sheds 
light  on  this  behavior.  We  find  that  (A R  monotonical- 
ly  increases  from  2.4  A  at  A  =  1 1  A  to  4.7  A  at  A  =  40  A, 
correlating  with  the  monotonic  increase  in  the  quality  of 
<r)  and  rm%.  The  fact  that  the  local  energy  is  relatively 
constant  at  large  r,  so  that  the  computed  energies  are 
only  weakly  dependent  on  the  accuracy  of  the  sampling 
in  this  domain,  readily  explains  the  difference  in  behavior 
of  the  energy  versus  the  coordinate  expectation  values. 

As  a  precursor  to  DMC,  we  also  performed  guided- 
walk  calculations  with  A  —  Dr  and  F=Vlni4'l2.  Note 
that  T  in  Eq.  13)  is  now  equivalent  to  Ga,  Eq.  (10),  if  the 
branching  is  omitted.  We  now  encountered  difficulty  in 
sampling  this  distribution  because  of  the  sharp  cutoff  at 
small  r.  In  computing  at  r=5Xl04,  10X104,  and 
20 X  104  hartree-1,  we  found  large  errors  in  the  energy, 
1%  for  the  first  two  time  steps  and  5%  for  the  last, 
despite  the  large  number  of  points  sampled, 
M  —  ( 3  — 6 )  X  107.  The  behavior  of  the  computed  values  of 
( r  )  and  r,mi  is,  however,  much  different.  Good  accuracy 
and  precision  are  obtained  for  these  quantities  at 
r—  10X  104  and  20X  104  hartree'1. 

The  reason  for  the  poor  estimates  of  the  energy  is  that 
walkers  are  either  ped  at  smaller  separations  or  else 
they  do  not  sample  these  domains.  This  trapping  is 
caused  by  the  guiding  force  F  being  excessively  large  at 
small  r,  where  it  is  proportional  to  r  "6,  giving  acceptance 
probabilities  practically  equal  to  zero.  This  in  turn  yields 
a  poor  representation  of  the  density  at  small  r  which  does 
affect  the  computed  energy  because  of  the  large  magni¬ 
tude  of  El  there.  This  is  much  less  significant  for  r  and 
rl  which  are  relatively  small  near  the  origin.  Since  the 
effect  of  F  on  the  acceptance  probability  is  removed  ex¬ 
ponentially  fast  as  r— 0,  the  small-r  domain  is  more  ac¬ 
curately  sampled  as  r  is  reduced.  However,  efficiency  is 
reduced  at  small  r  because  the  small  average  step  size 
gives  a  larger  degree  of  correlation  between  moves. 

We  now  consider  the  pure  helium  clusters.  These  sys- 
terns  possess  a  highly  repulsive  potential  as  does  lllle 
but  are  more  tightly  bound.  (For  example,  the  binding 
energy  of  He3  is  five  times  greater  than  that  of  I l-Me  ) 
This  causes  the  computation  of  expectation  values  by  an 
unguided  walk  to  be  less  difficult  than  for  H,Hc.  pus-urn- 
ably  because  sampling  a  slowly  decaying  distribution  at 
large  r  is  no  longer  required.  Despite  sampling  :."rr 
points  in  calculations  on  He,_3  7  [A/  =  ( 2  — 6) X  10*  w 
5X107  for  H.Hc],  for  each  cluster  we  observe  cvcclUul 
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agreement  between  computed  expectation  values  over  a 
range  of  A.  For  example,  with  2.6  <  A  <  13.2  A  for  He4, 
the  maximum  difference  in  ( r  >  (the  average  particle  sep¬ 
aration)  is  only  0.15%  and  that  in  the  energy  is  only 
0.07%. 

While  computed  averages  do  not  vary  significantly  as 
A  is  changed,  statistical  error  in  coordinate  expectation 
values,  ( r )  and  R =  R1)/])/)1'  ,  where 

K2  =  2;V-i<r,-Rc.m. >*].  decreases  as  A  (and  <A R))  in¬ 
creases.  In  going  from  the  smallest  to  the  largest  values 
of  A,  statistical  error  in  Rmi  and  (r )  is  continuously  re¬ 
duced  down  to  a  factor  of  2  or  more,  resulting  in  a  four¬ 
fold  or  greater  increase  in  computational  efficiency  for 
these  quantities.  For  the  energy,  small  values  of  A  result 
in  low  efficiency.  However,  once  A.  and  thereby  f  A R  ), 
is  of  reasonable  size,  statistical  error  in  the  energy  is  no 
longer  decreased  as  A  is  increased. 

Overall,  we  find  the  average  displacement  ( &R  )  to  be 
a  useful  measurement  of  the  sampling  efficiency  in  that 
larger  values  tend  to  give  smaller  statistical  errors,  most 
noticeably  for  R^  and  (r).  We  point  out  here  that  the 
sampling  required  to  obtain  <  A R  )  to  high  precision  is 
quite  small.  Therefore  finding  the  value  of  A  w  hich  yields 
the  greatest  average  displacement  may  be  accomplished 
with  very  little  computation.  Finally,  as  was  the  case 
with  H2He,  the  acceptance  ratios  corresponding  to  the 
largest  (A R  >  were  less  than  0.5,  i.e.,  0.38  for  He4  and 
0.40  for  He-,. 

We  now  turn  to  the  guided  walk,  which  encounters 
difficulty  when  sampling  at  small  r  for  H;He  and  does 
not,  therefore,  appear  useful  for  helium  clusters.  Howev¬ 
er,  since  the  DMC  walk  we  employ  consists  of  the  guided 
walk  described  above  with  branching,  evaluating  the 
practicality  of  this  guided  walk  is  important  in  ascertain¬ 
ing  the  feasibility  of  our  DMC  approach  for  obtaining  ex¬ 
act  results. 

As  discussed  in  Sec.  IV,  acceptance  probabilities  in  the 
guided  walk  globally  increase  as  the  time  step  is  reduced. 
If  a  time  step  can  be  found  which  is  small  enough  to  re¬ 
move  trapping  at  small  r,  so  that  convergence  in  sam¬ 
pling  ! ^ ! 2  can  be  obtained,  without  excessively  degrading 
sampling  efficiency,  DMC  may  be  practical  for  helium 
(and  other)  quantum  clusters.  Therefore  we  now  deter¬ 
mine  values  of  -  which  yield  high  accuracy  in  the  guided 
walk  for  He3_3. 

Table  II  presents  guided-walk  energies  and  statistical 
errors  (per  point  sampled)  for  He3_5  at  several  values  of  r. 
The  effect  of  trapping  is  immediately  evident  from  the 
data  in  Table  II.  At  the  larger  time  steps  the  energies  are 
of  poor  quality.  In  these  walks  we  have  observed  that 
atoms  which  are  too  close  together  do  not  move  during 
the  entire  course  of  a  run.  (Trapping  is  found  by  record¬ 
ing  the  number  of  accepted  moves  for  each  particle.)  As 
the  time  step  is  reduced,  particles  are  no  longer  trapped 
throughout  the  run.  The  result  is  a  noticeable  (and  for 
He,  a  dramatic)  improvement  in  the  energy.  Finally,  we 
see  that  at  sufficiently  small  time  steps,  guided-walk  ener¬ 
gies  are  in  excellent  agreement  with  those  computed  by 
the  unguided  approach. 

The  small-r  behavior  of  t he  sampling,  and  its  depen¬ 
dence  on  t,  is  illustrated  in  Fig.  1.  This  figure  compares 


TABLE  If.  Guided-walk  results  for  He,.,.  Time  steps  are 
given  in  hartree"1.  lengths  are  in  A.  and  energies  are  in  K. 


I0"4r 

<A  R) 

E/N 

olE/N) 

25.00 

2.14 

He, 

-0.0369(33) 

7.4 

12.50 

1.60 

-0.0367(27) 

6.0 

5.00 

1.09 

-0.0381(14) 

3.) 

2.50 

0.80 

-0.039  7(16) 

5.1 

1.25 

0.58 

-0.041  8(7) 

4.4 

0.75 

0.45 

-0.04141(12) 

1.1 

Unguided 

2.78 

-0.04147(7) 

0.14 

25.00 

1.86 

He4 

-0.1269(26) 

7.4 

15  00 

1.55 

-0.128  8(6) 

2.1 

10.00 

1.34 

-0.128  1(3) 

1.1 

5.00 

1.04 

-0.1329(23) 

92 

2.50 

0.78 

-0.136  3(19) 

10.7 

1.25 

0.57 

-0.136  1(7) 

4  4 

Unguided 

2.21 

-0.1356(1) 

0.2 

10.00 

1.17 

He, 

-0.34(10) 

200 

5.00 

0.96 

-0.34(10) 

200 

2.00 

0.69 

-0.248  7(12) 

3.4 

1.00 

0.51 

-0.25017(63) 

2.2 

0.50 

0.37 

-0.25048(51) 

2.5 

Unguided 

1.65 

-0.250  23(13) 

0.3 

r  (A) 


FIG.  1.  Convergence  of  pir)  at  small  r  in  the  guided  walk  for 

He,.  The  solid  line  ( - )  is  the  unguided  walk,  the  chain- 

dashed  line  ( - )  is  the  guided  walk  at  r=  50000  har*ree~'. 

the  long-dashed  line  ( - )  is  at  r  =  25  000  hartree'',  the 

slum-dashed  line  (-  -  -)  is  at  t—  12  500  hartree "  and  the  dotted 
line  (  -  ■  ■  )  is  at  r-7500  hartree"’.  Note  the  improving  agree¬ 
ment  between  uuguided  and  guided  as  the  time  step  for  the 
latter  is  reduced.  The  guided-walk  energies  follow  the  same 
trend,  see  Tjble  II. 
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the  probability  density  functions  p(r)  at  small  particle 
separations  (r)  corresponding  to  several  values  of  r,  with 
that  of  the  unguided  walk  for  He}.  Fluctuations  in  p  indi¬ 
cate  trapping  in  certain  regions  and  lack  of  penetration 
into  others.  These  fluctuations  decrease  as  r  is  reduced, 
resulting  in  convergence  to  the  unguided  probability  den¬ 
sity  function  and  agreement  in  the  computed  energies. 

It  is  important  to  point  out  that  the  error  in  the  energy 
at  the  larger  time  steps  is  not  systematic.  The  trapping  of 
points,  or  conversely  the  inability  to  sample  certain 
domains  at  small  particle  separation  (for  a  practical 
amount  of  sampling),  will  result  in  energies  either  too 
high  or  too  low,  depending  on  the  sign  of  the  local  energy 
at  small  r  and  whether  p  is  too  high  or  too  low.  For  the 
same  reasons,  statistical  error  may  be  artificially  high  or 
low  when  obtaining  an  accurate  representation  o{  p{r)  at 
small  r  is  problematic. 

For  purposes  of  comparing  efficiency  with  the  unguid¬ 
ed  walk,  the  statistical  error  in  the  energy  per  point  sam¬ 
pled,  <7,  is  presented  in  the  last  column  of  Table  II.  For 
M  points,  a  is  equal  to  the  statistical  error  in  the  average 
times  v  .Vf.  This  statistic  depends  on  the  sampling 
inefficiency  of  the  algorithm,  i.e.,  the  degree  to  which  suc¬ 
cessively  sampled  points  are  correlated,  as  well  as  the 
nonconstancy  of  EL,  which  is  determined  by  the  wave 
function.  Since  'V  is  the  same  for  the  computations  on  a 
given  cluster,  algorithm  efficiency  may  be  directly  com¬ 
pared.  Table  II  shows  that  while  agreement  with  unguid¬ 
ed  energies  is  obtained  by  the  guided  walk  at  sufficiently 
small  r,  the  unguided  walk  is  consistently  more  efficient, 
with  a  an  order  of  magnitude  smaller,  yielding  a  decrease 
in  computational  efficiency  of  two  orders  of  magnitude. 
Correspondingly,  for  the  guided  walk  at  the  smallest  time 
step,  the  average  step  size  {(SR  >)  is  four  to  five  times 
smaller  than  that  for  the  unguided  walk.  Nevertheless 
the  guided  walk  still  gives  good  precision  for  reasonable 
computational  effort.  If  this  is  also  the  case  with  the 
DMC  approach,  the  increase  in  accuracy  will  be  well 
worth  the  computational  effort. 

The  guided-walk  approach  we  have  employed  is 
inefficient  because  of  the  rapid  increase  in  the  guiding 
force  as  atoms  coalesce.  This  suggests  that  a  guided-walk 
approach  with  a  better  behaved  force  will  not  encounter 
the  difficulties  found  here.  We  have  investigated  two 
such  choices  in  a  few  selected  applications.  The  first, 
“weighted  unit  force,”  approach  simply  employs  cF,  in 
place  of  F,  =V,ln|'l/!i  when  moving  particle  i.  Thus  the 
direction  of  the  original  F  is  preserved  while  its  magni¬ 
tude  c  is  held  constant.  In  the  second,  a  “damped  force” 
approach,  y,F,  replaces  F,,  where 


*,=exp 


2  M'V/l 

j  i  *0 


(32) 


and  a0  in  /2  becomes  an  adjustable  parameter.  Once 
again,  the  direction  of  the  force  is  left  unchanged  but  now 
its  magnitude  is  most  greatly  affected  only  in  the  trapping 
regions,  i.e.,  !t,F,  |  — «0  when  V— .0. 

For  the  H2He  test  case  the  parameters  governing  the 
force,  c  and  o0,  arc  varied  together  with  r  to  obtain  max¬ 
imum  (  SR  >,  which  is  roughly  the  same  for  the  two  guid¬ 
ed  walks.  Trapping  is  not  observed  in  cither  case  and 


values  of  (SR)  are  obtained  which  are  a  factor  of  2 
greater  than  the  maximum  obtainable  in  the  unguided 
walk.  However,  the  statistical  error  in  the  energy  shows 
little  variation  among  the  unguided,  weighted,  and 
damped  approaches.  On  the  other  hand,  guided-walk 
statistical  errors  in  both  < r  >  and  rmi  are  about  25% 
lower  than  those  of  the  unguided  walk. 

For  the  second  test  case,  He!!2,  only  the  damped  force 
approach  was  compared  against  the  unguided  walk. 
While  exhibiting  no  trapping,  the  damped  force  walk 
yielded  no  increase  in  efficiency  over  the  unguided  walk. 
Therefore  we  conclude  that  the  simple  unguided  walk  is 
competitive  with  the  guided-walk  approaches  studied 
here. 


B.  DMC  Computations  on  helium  dusters 

For  the  DMC  walks,  inaccurate  sampling  at  small  r 
can  have  an  effect  significantly  greater  than  that  observed 
for  the  guided  walk.  If  points  are  temporarily  trapped  in 
a  region  where  the  local  energy  is  very  low,  as  we  have 
seen  can  easily  occur  at  small  r.  the  branching  factor  will 
be  very  large,  resulting  in  a  quickly  increasing  number  of 
walkers  at  small  r.  Although  the  particles  trapped  at  a 
small  separation  may  move  to  larger  r  after  a  short  period 
of  time,  the  continuous  generation  of  new  walkers  at  this 
point  will  yield  a  high  degree  of  oversampling  and  there¬ 
by  a  highly  biased  (too  low)  energy.  Therefore,  at  a  given 
time  step,  trapping  may  not  be  problematic  for  the  guid¬ 
ed  walk  but  nevertheless  gives  instability  in  the  DMC  ap¬ 
proach.  So  we  expect  that  smaller  time  steps  will  be  re¬ 
quired  to  obtain  convergence  in  the  DMC  walk  than  in 
the  guided  walk.  This  is  investigated  below. 

Given  the  instabilities  possible  in  DMC,  we  take  two 
steps  to  monitor  the  behavior  of  the  walk.  First,  we  al¬ 
low  the  ensemble  to  only  reach  twice  (or  four  times  for 
the  larger  clusters)  its  original  size.  If  the  maximum  en¬ 
semble  size  is  attained,  this  event  is  recorded  and  all 
weights  are  carried  so  that  the  copying  of  walkers  is  no 
longer  employed.  To  some  extent  this  step  controls  insta¬ 
bilities  arising  from  trapping  in  that  the  continuous  repli¬ 
cation  of  trapped  walkers  is  not  permitted.  Additionally, 
if  the  weight  of  any  walker  exceeds  a  given  value  (10),  the 
weight  is  recorded  in  the  output  file. 

Results  for  the  energy  and  its  growth  estimate  Ec  are 
presented  in  Table  III  and  plotted  in  Figs.  2(a)  —  2(c)  for 
Hej.5.  Block  length,  on  the  order  of  10s  — 106  hartree-1, 
was  varied  by  factors  of  2-4  resulting  in  no  significant 
change  in  the  computed  energies.  Also,  maximum  en¬ 
semble  sizes  (200!  and  excessive  weights  (10),  while  found 
at  the  larger  time  steps,  were  not  observed  at  the  smaller 
(last  two  or  three)  values  of  r.  Given  this,  and  the  statist¬ 
ical  agreement  of  the  energies  at  the  smaller  time  steps, 
the  DMC  energies  we  have  computed  are  deemed  to  be 
well  converged.  However,  we  do  see  that  much  smaller 
values  of  r  arc  indeed  required  than  was  the  case  for  the 
guided  walk;  approximate  comparisons  are  1000  versus 
12  500  hartree  1  for  He,,  1500  versus  25  000  hartree  1 
for  He4,  and  1500  versus  10000  hartree-1  for  He5  com¬ 
pare  Tables  II  and  III). 

Umrigar,  Runge,  and  Nightingale  have  recently  dc- 
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TABLE  III.  DMC  energy  vs  time  step  for  He,.,.  The  average  of  the  local  energy  is  E  and  the 
growth  estimate  is  Ec.  Time  steps  are  given  in  hartree*1,  lengths  are  in  A,  and  energies  are  in  K.  AR 
denotes  the  acceptance  ratio. _  _ 


o 

1 

w 

AR 

<ER) 

E/N 

Ec/S 

4.0 

0.9963 

0.62 

He, 

-0.048  6(14) 

0.044  7(7) 

3.0 

0.9975 

0.54 

-0.046  2(7) 

-0.044  28(35) 

2.0 

0.9986 

0.44 

-0.044  33(16) 

-0.044  06(17) 

1.0 

0.9995 

0.31 

-0.04409(18) 

-0.04398(17) 

0.5 

0.9998 

0.22 

-0.04428(20) 

-0.043  93(23) 

5.0 

0.9908 

0.69 

He, 

-0.1750(90) 

-0.145  5(7) 

2.5 

0.9965 

0.49 

-0.1514(38) 

-0.145  4(6) 

2.0 

0.9975 

0.44 

-0.148  F  16) 

-0.145  7(7) 

1.5 

0.9983 

0.38 

-0.144  5(3) 

-0.144  8(4) 

1.0 

0.9991 

0.31 

-0.144  5(3) 

-0.144  4(3) 

2.0 

0.9963 

0.41 

He, 

-0.282(13) 

-0.268  2(8) 

1.5 

0.9975 

0.38 

-0.268  3(5) 

-0.267  5(4) 

1.0 

0.9986 

0.31 

-0.267  3(11) 

-0.2674(10) 

scribed  a  DMC  approach  which  controls  the  magnitude 
of  F  near  its  singularities  and  yields  a  better  approxima¬ 
tion  to  the  Green’s  function  in  this  domain  [41].  It  will 
be  of  interest  to  see  if  this  method  will  reduce  time-step 
bias,  and  thereby  increase  efficiency  by  allowing  greater 
values  of  r,  in  DMC  computations  on  helium  clusters. 

Neglecting  the  different  time-step  scales,  the  behavior 
of  the  DMC  energies  is  very  similar  to  those  of  the  guid¬ 
ed  walk.  In  essence,  both  walks  are  affected  by  trapping, 
which  is  magnified  by  the  branching  in  the  DMC  walk. 
As  was  also  the  case  for  the  guided  walk,  the  coordinate 
expectation  values,  (r)  and  Rrm,  were  not  visibly 
affected  by  r.  This  is  not  surprising  as  these  quantities 
are  only  weakly  influenced  by  sampling  accuracy  at  small 
particle  separations.  This  lack  of  time-step  bias  is  also 
seen  for  the  particle  separation  density  functions  p(r)  (ex¬ 
cept  of  course  at  small  r),  as  well  as  for  the  density 
profiles. 

Computational  cost  in  obtaining  converged  DMC  ener¬ 
gies  was  large  but  not  excessive.  While  He,  presented  an 
especially  difficult  case  requiring  five  Cray  X-MP14  CPU 
hours  for  all  (not  each)  of  the  time  steps,  He4  and  He5 
took  only  one  and  two  hours,  respectively.  As  for  H;He 
at  the  VMC  level  of  theory,  the  smallest  and  most  diffuse 
cluster  gave  the  greatest  difficulty  in  obtaining  the  accu¬ 
racy  and  precision  desired. 

Having  successfully  obtained  converged  energies  for 
the  smaller  clusters,  it  was  of  interest  to  see  if  this  could 
also  be  achieved  for  the  larger  clusters.  We  found  for 
Hev,  N  =  7,  14,  20,  and  1 12,  DMC  energies  converged  at 
about  t—  500-  1500  hartree”1.  Even  at  the  largest  time 
step  of  2000  hartree"1,  the  energies  were  very  close  to 
those  computed  at  the  smallest  values  of  t.  Overall,  the 
dependence  of  F  /N  and  Ec/S  versus  the  time  step  mim¬ 
ics  that  of  the  smaller  clusters.  At  “larger”  r  (  =2000 
hartree  ’),  bias  is  noticeable  and  large  weights  and  fluc¬ 
tuations  in  ensemble  size  occur.  At  ".smaller"  r  (  ~  1000 


hartree-1),  stability  in  the  DMC  energies,  weights,  and 
ensemble  sizes  is  obtained.  Therefore,  while  the  required 
time  step  for  unbiased  energies  is  greatly  reduced  for 
DMC  versus  the  guided  walk,  this  is  not  the  case  for 
larger  clusters  versus  smaller  clusters.  Unbiased  energies 
are  obtained  at  r-  1000- 1500  hartree-1  for  3-5  atoms 
3nd  t  =  50G-1500  hartree-1  for  7,  14,  20,  and  112  atoms. 
(Here,  absence  of  bias,  i.e.,  it  being  masked  by  statistical 
error,  is  relative  rather  than  absolute.  For  example,  a 
bias  of  0.005  K  is  very  large  for  He},  1 1%  of  E  AV,  but  on 
the  order  of  the  statistical  error  for  Hn;,  0.1%  of  E/S.) 
This  is  explained  by  the  fact  that  convergence  in  DMC  is 
most  affected  by  trapping  combined  with  large  fluctua¬ 
tions  in  the  local  energy  at  small  r  —  effects  which  are 
governed  mostly  by  the  wave-function  form  rather  than 
by  cluster  size.  Total  computational  cost  was  roughly  2, 
4,  5,  and  17  Cray  X-MP14  CPU  hours  for  N  =  7,  14,' 20, 
and  112,  respectively. 

The  DMC  results  are  summarized  and  compared  to 
GFMC  [10,11]  and  other  recent  DMC  results  [18]  (He;0 
and  Hen2)  in  Table  IV.  The  results  we  compare  against 
were  obtained  with  the  HFDHE2  potential  [32]  which 
predates  the  most  recent.  HFD-B(HE),  potential  [31] 
used  here.  The  two  potentials  possess  practically  identi¬ 
cal  functional  forms  but  with  different  sets  of  parameters. 
Perhaps  most  significant  is  the  f  .3%  increase  in  well 
depth.  For  the  unit  radius,  r0  =  v/5/3Rmu/.Vl/\  we  em¬ 
ploy  a  "second-order”  approximation  of  exact  expecta¬ 
tion  values  [7]  defined  as 

<  A  >Is2<'PM!<ft0)/<'Pl<i0)-<'P[/Ii4'>  ,  (33) 

with  A  =R2.  Writing  'l'  =  <5u  +  5  shows  that  (  A  >,  differs 
from  the  exact  value,  (<5„| .-I ]<50),  by  integrals  involving 
6'.  (This  approximation  is  useful  when  A  does  not  com¬ 
mute  with  the  Hamiltonian.  Methods  for  computing 
( i30|  .1 ; c/>0 )  have  been  described  elsewhere  [42,43].) 
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The  DMC  energies  we  compute  with  the  HFD-B(HE) 
potential  are  significantly  below  the  GFMC  energies  re¬ 
sulting  from  the  previous,  HFDHE2,  potential.  The  rela¬ 
tive  discrepancies  tend  to  decrease  with  increasing  cluster 
size.  For  He3  our  energy  is  13%  lower  than  the  GFMC 
value  while  for  He,*  and  He,0  the  differences  are  only 
3.3%  and  3.7%,  respectively.  However,  as  seen  in  Table 
IV,  for  He3  ana  He2o  cur  DMC  energies  are  in  excellent 
agreement  with  GFMC  when  the  same  potential 
sHFDHE2)  is  employed.  In  contrast,  our  HFDHE2  ener¬ 
gy  does  not  agree  with  GFMC  for  Helu.  The  new  poten¬ 
tial  lowers  our  Hell2  DMC  energy  by  3.2%.  Such  sensi¬ 
tivity  to  small  changes  in  the  potential  has  been  observed 
previously.  Kalos  et  al.,  employing  the  HFDHE2  poten¬ 


tial  in  their  study  of  liquid  4He,  obtained  a  6%  decrease 
in  the  energy  for  a  1.9%  increase  in  the  well  depth  (44). 

Upon  considering  He,t2  and  comparing  with  the  DMC 
energies  computed  by  Chin  and  Krotscheck  [18]  (CK), 
employing  a  DMC  algorithm  different  than  our  own  [45], 
discrepancies  not  attributable  to  the  potential  arise.  For 
He20,  we  see  in  Table  IV  that  CK’s  energy  is  2%  below 
both  our  DMC  and  the  GFMC  values,  when  ail  three  cal¬ 
culations  employ  the  same  (HFDHE2)  potential.  For 
Hen2,  further  disagreement  occurs  with  the  same  poten¬ 
tial  as  our  and  CK's  energy  lie  roughly  2%  and  3%,  re¬ 
spectively,  below  GFMC.  This  is  not  readily  explained 
by  statistical  error  which,  for  both  GFMC  and  DMC,  is 
an  order  of  magnitude  smaller  than  these  differences.  It 


FIG.  2.  Local  (circles)  and  growth  (squares)  energy  estimates  of  the  exact  energy  vs  r  for  (a)  He,,  lb)  He.,,  and  le)  He,. 
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tion  has  equilibrated.  This  was  done  for  He3,  reducing  same  potential.  The  difference  between  our  VMC  and 
the  discrepancy  between  the  two  energy  estimates  from  DMC  values  is  only  0.6%,  implying  an  error  of  much  less 

0.9%  to  0.5%,  thereby  giving  near  statistical  agreement.  than  0.6%  in  the  second-order  estimate.  Statistical  error 

For  the  unit  radii  r0,  we  first  compare  our  DMC  and  in  the  second-order  values,  0.003  A  or  less,  cannot  be  the 

the  GFMC  values  for  He3  and  HeJ0  computed  using  the  cause  of  the  2%  discrepancy  between  second-order  DMC 

same  potential  (Table  IV).  While  agreement  is  observed  and  GFMC. 

for  He:0,  a  large  difference  of  8%  is  seen  for  Hes.  This  The  unit  radii  computed  by  CK  differ  slightly  from  our 
difference  is  well  beyond  our  statistical  error  of  1%.  (Ex-  own,  1.5%  below  fer  Hc2q  and  0.8%  above  for  He,|2. 

cept  for  He3,  statistical  error  in  DMC  values  of  rQ  is  well  These  differences  are  most  likely  caused  by  the  differences 

under  1%,  i.e.,  0.01-0.02  A.)  Furthermore,  this  in  the  DMC  solutions  obtained,  as  reflected  by  the  ener- 

disagreement  with  GFMC  does  not  appear  to  be  caused  gies,  rather  than  by  the  second-order  approximation, 

by  error  in  our  second-order  estimate  of  r0.  The  error  in  This  is  supported  by  the  fact  that  CK  also  obtain  very 

the  second-order  estimate  of  r0,  z=(8\R2\5),  should  be  good  agreement  between  their  VMC  and  second-order 

less  than  the  difference  between  the  VMC  and  DMC  unil  radii. 

values.  R-”> {'),  which  is  only  2%.  For  the  largest  Finally,  in  comparing  results  using  the  HFD-B(HE) 

cluster,  N  =  112,  disagreement  resurfaces  when  compar-  potential  for  He)U,  note  that  the  difference  between 

ing  our  value  of  r0  against  GFMC  obtained  with  the  VMC  and  DMC  is  only  0.6%.  Therefore,  in  this  case,  we 


R  (A)  R  (A) 


HG.  3,  VMC  ;imi  second -order  density  profiles  for  (a)  He,,  r=500  hartree  (W  HC|4,  r  =  5 00  hartree  '.  (c)  He.-#,  r  “750  h.ir- 
tree  \  ,uid  <d;  He,,.,  r  -  HXX)  hartree  The  dotted  line  is  the  VMC  and  the  solid  line  is  the  second -order  approximation  to  the  ex¬ 
act  ;!iq.  (33iJ.  The  straight  line  in  id)  represents  the  experimental  liquid-helium  density  ol  0  021  85  A  ~J. 
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also  expect  our  computed  value  of  r0  to  be  accurate  to 
the  number  of  figures  shown.  The' change  in  r0  on  com¬ 
paring  the  old-  versus  the  new-potential  values  (2.396 
versus  2.390  k)  is  only  0.3%,  significantly  less  than  the 
1%  difference  we  find  for  He20.  In  all  cases  where  com¬ 
parison  is  made,  r0  is  reduced  when  employing  the  poten¬ 
tial  with  the  deeper  well,  as  expected.  This  effect  de¬ 
creases  by  an  order  of  magnitude,  3.8%  to  0.3%,  in  pass¬ 
ing  from  three  to  1 12  atoms  and  contrasts  with  the 
change  in  the  energy,  13.2%  to  3.2%. 

Of  greater  interest  is  the  reliability  of  our  reported  esti¬ 
mates  of  r0  obtained  with  the  most  up-to-date  potential. 
We  have  found  that  the  differences  between  VMC  and 
DMC  reach  a  maximum  of  only  3%  (He4  and  He7). 


Therefore  we  expect  an  accuracy  of  at  least  3%  in  our  es¬ 
timates  of  r0.  However,  the  accuracy  is  probably  much 
better;  for  He3  0.5%  agreement  is  obtained  from  wave 
functions  with  noticeably  different  VMC  values  of  r0, 
5.50  and  5.74  A,  i.e.,  above  and  below  the  estimated  exact 
value. 

Figures  3(a)- 3(c)  resent  VMC  and  second-order 
(r=5C0- 1000  hartree"1)  [see  Eq.  OB')]  density  profiles 
for  He^,  HeM,  and  He20.  The  bin  sizes  are  very  small, 
0.026  A  for  He3  and  0.019  A  for  Heu  and  He^.  This  al¬ 
lows  for  very  fine  detail  in  p,  and  although  there  is 
enhanced  statistical  error  near  the  origin,  it  remains 
below  10%  down  to  about  0.2-1  A.  We  see  for  these 
clusters,  as  well  as  for  the  other  clusters  not  shown,  that 


FIG.  4.  VMC  and  second-order  probability  density  functions  of  particle  separation  for  (a)  He,,r=500  hartree'1,  (b)  He,4.  r  =  500 
hartree  " 1 ,  and  (e)  He20,  r=750  hartree'1.  The  dotted  line  is  the  VMC  and  the  solid  line  is  the  second-order  approximation  to  the  ex¬ 
act. 
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the  qualitative  behavior  is  unchanged  between  the  VMC 
and  second-order  profiles.  The  sharp  increase  in  p  near 
the  cluster  center  for  He3,  which  we  first  observed  at  the 
VMC  level  of  accuracy  and  which  arises  from  a 
significant  contribution  of  near-collinear  configurations 
to  the  density  [28],  remains  in  DMC.  For  14-  and  20- 
atom  clusters,  very  little  structure  is  evident.  The  Heu 
density  rises  slightly  near  the  origin  while  that  of  He20 
reaches  a  constant  value  at  about  2.5  k  and  then  fluctu¬ 
ates  about  0.019  A  3.  This  is  in  good  agreement  with 
GFMC  [10]  and  the  oscillations  seen  by  CK  are  not  ob¬ 
served  here. 

VMC  and  second-order  (r=500  hartree-1)  density 
profiles  for  Hem  are  computed  with  a  bin  size  of  0.13  A 
and  are  presented  in  Fig.  3(d).  The  experimental  liquid- 
helium  density  of  0.021  85  A  3  is  shown  for  comparison 
as  a  solid  line.  Statistical  error  is  10%  for  the  points 
nearest  the  origin,  then  decreases  rapidly  at  greater  dis¬ 
tances,  and  finally  begins  to  rise  near  the  cluster  edge, 
reaching  about  10%  at  1 1  A.  Unlike  any  of  the  smaller 
clusters,  structure  in  the  density  profile  now  appears  to  be 
present.  However,  further  analysis  indicates  that  the 
fluctuations  at  R  <  5  A  are  statistical.  Only  in  a  very 
small  region,  2.10-2.35  A,  is  statistical  error  (4.5%) 
significantly  below  the  deviation  from  the  liquid  density 
(6% -10%).  Thus,  out  to  4.8  A,  the  density  is  very  nearly 
constant  to  within  reasonable  statistical  error  (under  5% 
for  R  >  1.5  A).  However,  a  shoulder  is  present  at  6.2  A 
where  statistical  error  is.  very  small,  1.5%.  Further  out, 
another  shoulder  at  9.6  A  is  barely  discernible.  This  is  in 
good  agreement  with  CK  <b;n  size  is  0.24  A),  who  also 
observe  shoulders  in  these  regions.  We  have  also  com¬ 
puted  a  second-order  density  profile  for  HeU2  employing 
the  HFDHE2  potential.  The  results  are  very  similar  to 
those  we  obtained  above. 

In  summary,  our  second-order  densities  for  N  —  14,  20, 
and  112  rise  up  to  their  central  values  at  some  distance 
from  the  origin.  The  value  of  this  central  density  and  the 
distance  to  which  it  extends  increases  with  cluster  size. 
The  Heu  central  density  is  clearly  below  that  of  liquid 
helium  and  p  begins  to  drop  off  at  about  1  A.  That  of 
He,  |2  is  in  good  agreement  with  the  liquid-He  density 
and  extends  out  to  about  4.5  A  while  the  He;0  case  is  in¬ 
termediate  between  Heu  and  He,u.  These  conclusions 
differ  from  those  of  CK,  who  obtain  oscillations  in  their 
density  profiles  as  the  origin  is  approached  for  both  He,0 
and  He,l2.  We  see  no  such  oscillations  for  He:o  while 
those  of  Hen2  are  mostly  removed  upon  considering  the 
small  (2%-5%)  statistical  error.  However,  we  do  see 
shoulders  in  the  Hell2  density  at  larger  R. 

For  all  clusters,  the  quantitative  differences  between 
the  VMC  and  second-order  density  profiles  are  not  large. 
The  small  changes  we  observe  in  passing  from  VMC  to 
second  order  lead  us  to  believe  our  estimate  of  exact  den¬ 
sity  profiles  by  second-order  profiles  is  well  jusiified  (see 
discussion  concerning  the  unit  radii).  Comparison 
against  density  profiles  computed  by  VMC  [26]  or 
GFMC  [10]  shown  little  variation  in  the  central  density, 
despite  the  different  potentials  and  wave  functions  em¬ 
ployed.  A  major  difference  does  arise  for  He,  which  has 


FIG.  $.  VMC  ( - i  and  DMC  plots  cf  p(r)  at  short  range 

for  He14  at  r=2000  ( - 1,  1000  ( - I,  and  500  I-  -  -i  har¬ 

tree 


a  sharp  increase  in  p  at  about  2  A  [28]  which  is  not  evi¬ 
dent  in  GFMC  [10]. 

In  order  to  gain  further  insight  into  cluster  structure 
and  its  changes  as  accuracy  is  improved  from  VMC  to 
the  second-order  approximation,  we  compute  particle 
separation  probability  density  functions  p(r).  Figures 
4'a)-4(cl  present  VMC  and  second-order  plots  o (plr)  for 
He7,  He, 4,  and  He20.  For  He3_5,  qualitative  distinctions 
are  not  discernible  between  the  VMC  and  second-order 
density  functions,  just  as  for  He,.  For  those  clusters 
which  show  structure  in  p  with  VMC,  He14  and  He-,,, 
qualitative  differences  between  VMC  and  second-order 
densities  are  now  apparent.  We  see  that  the  implied 
shoulder  in  p  found  by  VMC  becomes  much  more  pro¬ 
nounced  as  we  progress  to  the  second-order  level  of  accu¬ 
racy.  It  would  be  interesting  to  see  if  this  onset  of  “shell’’ 
structure  is  progressive  or  abrupt  as  cluster  size  is  in¬ 
creased  from  7  to  14  atoms.  We  point  out  once  more  that 
both  p  and  p  are  essentially  independent  of  the  time  step 
employed  (for  the  range  of  r  we  have  considered).  The 
exception  of  course  is  for  p  at  small  r.  An  example  of  this 
is  given  in  Fig.  5  which  presents  the  VMC  and  several 
second-order  estimates  in  the  region  p(r)/pm„  <5X10  1 
for  Heu.  Note  that  by  r—  1000  hartree'1  convergence  is 
quite  good,  as  is  the  case  for  the  energy.  (See  Table  IV 
and  associated  discussion.) 

C.  New  two-body  form  and  T} 

We  have  employed  our  new  two-body  form,  t\  [see  Eq. 
(24)],  in  unguided-walk  computations  on  Hev,  .Y  =  J,  7, 
and  20.  For  each  cluster  the  same  value  of  A  is  employed 
for  / 2  and  t\.  We  find  that  f,  yields  the  best  energy  for 
He3,  5%  below  that  of  f2.  However,  for  the  larger  dus¬ 
ters,  small  but  significant  (i.e.,  well  beyond  statistie.il  er¬ 
ror)  reductions  are  obtained  with  i.e.,  2%  for  He,  and 
1%  for  He:o.  For  He3,  the  most  diffuse  cluster,  the  re- 


4096 


R.  N.  BARNETT  AND  K.  B.  WHALEY 


47 


duced  flexibility  of  t\  in  descnbing  jong-range  behavior  ts 
significant.  In  addition  to  yielding  a  poorer  energy,  i'2 
also  gives  an  unrealistically  low  unit  radius.  For  the 
more  tightly  bound  clusters,  for  which  short-range  in¬ 
teractions  are  more  important,  an  improvement  with  is 
obtained.  Also,  /2  yields  a  30%  reduction  in  statistical 
error  in  the  energy,  even  for  He,,  presumably  by  better 
descnbing  the  short-range  behavior  where  the  local  ener¬ 
gy  possesses  its  greatest  fluctuations.  Consequently,  we 
expect  that  further  accuracy  may  be  obtainable  by  in¬ 
creasing  the  long-range  flexibility  of  f  2. 

In  addition  to  seeking  better  two-body  wave  functions, 
accuracy  can  also  be  increased  by  including  three-body 
f3nd  higher)  effects,  as  discussed  in  Sec.  IV  A.  We  treat 
here  the  larger  clusters,  whose  wave  functions  should 
possess  the  greatest  need  for  three-body  correlations.  For 
He,4  and  He,0,  the  initial  parameters  in  7\  (/,  here)  are 
obtained  from  a  conjugate-gradient  optimization.  For 
Hen;  mi*'3!  parameters  are  those  for  He,0  and  are  then 
varied  by  hand.  As  seen  in  Table  I,  optimization  resulted 
from  changing  only  the  long-range  parameters.  Upon  ad¬ 
dition  of  Tx,  parameters  are  varied  by  hand  in  a  set  of 
short  Monte  Carlo  runs.  At  this  stage,  only  variation  of 
the  7~,  parameters  and  of  the  long-range  parameters  in 
T;  was  found  to  be  fruitful.  Despite  the  approximate  lev¬ 
el  of  optimization,  a  significant  reduction  in  the  energy  is 
observed  in  all  cases.  The  optimized  T}  term  yielded 
about  a  6%  improvement  in  the  energy  for  Heu  and 
He;o.  As  expected,  the  Hen:  energy  is  reduced  by  a 
greater  amount  (9%)  than  for  the  smaller  clusters.  The 
final  result  is  that  the  VMC  energies  of  the  .V  =  14,  20, 
and  112  clusters  are  quite  good;  96. 1  ^  of  the  computed 
exact  value  is  obtained  for  Hel4  and  this  decreases  only 
by  1.5 %  upon  going  to  Hen2. 

Wave-function  quality  is  also  improved  in  other 
respects.  The  data  listed  in  the  middle  of  Table  V  show 
that  the  relative  statistical  error  in  the  energy  decreases  in 
every  case.  The  increased  computation  when  T,  is  in  the 
wave  function  ranges  from  82%  for  Heu  to  88%  for 
Hen2.  Allowing  for  this,  the  efficiency  (the  inverse  prod¬ 
uct  of  the  variance  and  the  computation  time)  in  comput¬ 
ing  the  energy  is  increased  by  a  factor  of  2  for  Heu  and 
Hen2.  Interestingly,  however,  efficiency  is  decreased  by 
19%  for  He;o.  This  contrast  with  the  Heu  and  He,t2 
cases  may  arise  from  incomplete  (hand)  parameter  optim¬ 
ization  or  from  use  of  the  energy,  rather  than  variance, 
minimization  criterion.  It  could  also  be  magnified  by  the 
generally  large  uncertainties  in  computed  statistical  er¬ 
ror. 

TABLE  V.  Accuracies  and  precisions  for  two-  and  ihrce- 
hody  wave  functions. 

Accuracy  in  E/X  </(  £/.V)/(£Y,V)  Accuracy  in  rn 
(%)  (%) 


,v 

t2 

r.  +  r, 

Ti 

T;  +  Tx 

r. 

Tj  +  T 

14 

90.5 

96.1 

1  00 

0  52 

3.2 

1.2 

20 

89  5 

95.3 

0.93 

0.76 

4  7 

2.8  , 

112 

85  7 

94.6 

0  23 

0  II 

4.9 

1.2 

In  the  last  two  columns  of  Table  V,  percent  differences 
between  VMC  and  second-order  unit  radii  r0  are  listed 
for  the  7\  and  7"2  +  r3  wave  functions  of  Heu,  He;o,  and 
He,,,.  We  expect  these  differences  to  be  good  estimates 
of  the  deviation  rrom  exact  values,  giver,  our  confidence 
in  the  accuracy  of  our  second-order  estimates  of  r0,  see 
Sec.  V  B.  In  each  case  agreement  with  second-order  esti¬ 
mates  of  r0  is  noticeably  enhanced  upon  addition  of  the 
three-body  correlation  functions.  The  result  is  that  very 
good  agreement  (  =1%)  with  estimated  exact  unit  radii  is 
obtained  at  the  VMC  level  of  theory  with  the  three-body 
wave  functions.  The  exception  is  again  He:0  with  a  VMC 
value  of  rG  differing  by  3%  from  the  second-order  esti¬ 
mate. 

VI.  CONCLUSIONS 

We  have  studied  wave-function  forms  and  Monte  Car¬ 
lo  integration  techniques  for  H,He  and  Hev, 
.V  =3-20,112.  As  a  very  diffuse  system  with  a  highly 
repulsive  potential,  H-He  presented  special  difficulties  for 
the  VMC  computations.  While  the  VMC  approaches 
used  here  are  without  bias,  errors  in  computed  quantities 
can  arise  for  a  finite  (yet  large)  amount  of  sampling  if 
efficiency  is  sufficiently  poor.  Therefore  only  for  very 
large  attempted  displacement  sizes  A  is  good  accuracy 
obtained  in  the  unguided  walk.  Furthermore,  these 
values  of  A  correspond  to  acceptance  ratios  of  0.35-0.40, 
well  below  the  often  assumed  optimum  of  0.5. 

Despue  the  increased  dimensionality  in  comparison  to 
H;He  (3.V  versus  3).  the  less  diffuse  helium  clusters  are 
much  more  amenable  to  Monte  Carlo  integration.  For 
the  unguided  walk,  consistency  in  computed  expectation 
values  is  obtained  over  a  wide  range  of  A.  However,  we 
do  find  that  statistical  error  in  bolh  (r)  and  in  R^  is 
noticeably  lower  at  large  A.  Once  again,  these  values  of 
A  corresponded  to  the  largest  average  displacement  of  an 
aitempted  move.  (  A R  ),  rather  than  to  an  acceptance  ra¬ 
tio  of  0.5. 

Directly  including  importance  sampling  by  using  a 
walk  guided  by  F  =  V  'Fi'  yielded  inaccurate  H:He  ener¬ 
gies.  Although  most  of  configuration  space  was  well  sam¬ 
pled.  as  indicated  by  accurate  values  of  (r)  and  rmt, 
large  values  of  the  force  F  caused  trapping  at  small  r. 
The  resulting  inaccuracy  in  the  density  gave  rise  to 
significant  errors  in  the  computed  energy  because,  al¬ 
though  the  wave  function  is  small,  the  large  magnitude  of 
ihe  local  energy  at  small  separations  requires  accurate 
sampling  in  this  region.  However,  trapping  is  readily  cir¬ 
cumvented  by  employing  a  better  behaved  force.  Two 
such  approaches,  which  siill  direct  moves  toward  a  local 
maximum  in  I'Vl2,  were  applied  to  HTIe  and  Hen2. 
Compared  to  an  unguided  walk,  efficiency  was  increased 
slightly  for  H,He  but  not  for  Hen2. 

As  seen  for  H:He.  errors  in  the  energy  arise  in  a  walk 
guided  by  F  =  V  In  4,|‘  due  to  poor  convergence  at  small 
r  (which  is  practically  infinitely  long  unless  r  is  very 
small)  In  light  of  the  desire  to  compute  exact  energies  by 
DMC.  however,  it  was  of  interest  to  determine  (lie 
domain  of  r  at  which  high  accuracy  could  be  obtained  by 
ibis  guided  walk.  We  find  ihat,  though  initially  quite 
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poor  at  large  r,  accurate  energies  are  obtained  at  smaller 
values.  “Convergence"  in  the  guided-walk  energies  is 
directly  correlated,  as  expected,  with  improving  accuracy 
in  p{r)  at  small  r  as  the  time  step  is  reduced.  This  was 
also  the  case  in  the  DMC  calculations,  demonstrating 
that  accuracy  in  the  (DMC  or  VMC)  energy  is  critically 
dependent  on  the  sampling  at  small  r.  In  light  of  this  dis¬ 
cussion,  we  point  out  that  a  DMC  approach  similar  to 
that  employed  here  but  which  bounds  the  magnitude  of  F 
near  4'  =  0  has  been  described  recently  [41]  and  may  be  of 
use  for  the  systems  studied  here. 

Although  small  time  steps  are  required,  well-converged 
DMC  energies  have  been  obtained  for  Hev. 
.V  =3-5,7, 14,20, 1 12.  The  steps  taken  to  ascertain  con¬ 
vergence,  variation  of  the  time  step  and  block  size,  com¬ 
parison  of  the  local  and  growth  energy  estimates,  and 
convergence  of  p  at  small  r,  all  support  the  reliability  of 
the  computed  energies.  As  we  have  already  found  with 
VMC  for  He3  and  He4  [28],  energies  well  below  those  of 
OFMC  [10,11]  are  obtained.  For  He3  and  He;o,  these 
discrepancies  are  clearly  caused  by  the  use  of  different 
potentials.  It  is  reasonable  to  conclude  that  this  is  also 
the  case  for  all  the  3-20  atom  clusters  studied  here.  The 
agreemetu  we  obtain  with  GFMC  for  an  identical  poten¬ 
tial  leads  us  to  believe  in  the  reliability  of  our  DMC  ap¬ 
proach  and  the  exactness  of  our  computed  energies  with 
the  most  up-to-date  potential.  However,  energies  below 
GFMC  (but  with  the  same  potential)  have  been  obtained 
by  CK  using  a  modified  DMC  approach  for  He  v,  .V  =20, 
40,  70,  and  112  [18],  and  by  us  for  N— 112.  While  the 
discrepancy  between  DMC  and  GFMC  is  smaller  in  our 
calculation  than  in  CK’s,  our  energy  is  still  significantly 
lower  than  the  GFMC  result.  Further  disagreement 
occurs  when  comparing  our  DMC  energies  with  those  of 
CK.  For  both  ,V=20  and  112,  we  compute  slightly 
higher  energies. 

We  see  a  sizable  lowering  of  the  energy  below  that  ob¬ 
tained  from  the  earlier,  HFDHE2,  potential  when  em¬ 
ploying  the  most  recent,  HFD-B(HE),  He-He  interaction 
potential,  13%  for  three  atoms  and  3.2%  for  112.  This 
reflects  primarily  the  lower  well  depth  of  the  newer  po¬ 
tential.  Decreases  in  the  unit  radius  are  also  observed, 


3.8%  for  He3  to  0.3%  for  Hen2.  Finally,  as  do  CK,  we 
see  fluctuations  in  the  Hel)2  density  profile  which  have 
not  previously  been  observed  at  either  the  VMC  or 
GFMC  level  of  theory.  However,  our  fluctuations  at 
small  R  are  beneath  statistical  error  as  differences  from 
the  liquid-helium  density  are  generally  less  than  one  stan¬ 
dard  deviation  in  this  region,  R  <  5  A. 

In  an  effort  to  improve  accuracy  at  the  two-body  level, 
we  have  studied  an  entirely  new  form  describing  these 
etfects.  This  form  gives  added  emphasis  at  small  r  and 
contains  a  factor  which  mimics  the  potential  in  this 
domain.  Optimized  wave  functions  for  He7  and  He^ 
gave  slightly  improved  energies,  despite  the  reduced  flexi¬ 
bility  at  large  r.  For  the  more  diffuse  He3,  the  older  form 
gave  a  lower  energy.  In  addition,  statistical  error  in  the 
energy  was  reduced  by  about  a  third.  It  is  expected  that 
a  better  description  of  the  long-range  behavior  will  yield 
further  improvements. 

In  order  to  investigate  the  accuracy  obtainable  by 
current  VMC  approaches,  a  three-body  factor  was  added 
to  the  14-,  20-,  and  112-atom  two-body  wave  functions. 
Substantial  improvement  in  the  energy  is  obtained,  and 
for  He, 4  and  Hen2  an  increased  efficiency  in  computing 
this  quantity  also  results,  despite  the  greater  complexity 
of  the  wave  function.  More  sophisticated  optimization 
algorithms  for  the  parameters  in  T]  may  yield  a  further 
lowering  of  the  VMC  energy.  It  also  remains  to  be  seen 
whether  t':,  combined  with  T3,  will  yield  better  agree¬ 
ment  with  exact  energies,  and  whether  the  use  of  such 
complex  wave  functions  will  be  advantageous  for  DMC. 
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